Author: [log in to unmask]
Date: Wed Sep 2 17:41:06 2015
New Revision: 3501
Log:
compute covariance matrix, maybe
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/HpsGblRefitter.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
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 Wed Sep 2 17:41:06 2015
@@ -160,8 +160,8 @@
// Use the truth helix as the initial track for GBL?
//htf = htfTruth;
// Get perigee parameters to curvilinear frame
- PerigeeParams perPar = new PerigeeParams(htf);
- PerigeeParams perParTruth = new PerigeeParams(htfTruth);
+ PerigeeParams perPar = new PerigeeParams(htf, _B.z());
+ PerigeeParams perParTruth = new PerigeeParams(htfTruth, _B.z());
if (textFile != null) {
textFile.printPerTrackParam(perPar);
textFile.printPerTrackParamTruth(perParTruth);
@@ -171,9 +171,9 @@
gtd.setPerigeeTrackParameters(perPar);
// Get curvilinear parameters
- ClParams clPar = new ClParams(htf);
- ClParams clParTruth = new ClParams(htfTruth);
if (textFile != null) {
+ ClParams clPar = new ClParams(htf, _B.z());
+ ClParams clParTruth = new ClParams(htfTruth, _B.z());
textFile.printClTrackParam(clPar);
textFile.printClTrackParamTruth(clParTruth);
@@ -812,27 +812,27 @@
return Math.sqrt(Math.pow(E1 + E2, 2) - VecOp.add(p1vec, p2vec).magnitudeSquared());
}
- private BasicMatrix getPerParVector(HelicalTrackFit htf) {
+ private static BasicMatrix getPerParVector(HelicalTrackFit htf, double B) {
BasicMatrix perPar = new BasicMatrix(1, 5);
if (htf != null) {
- double kappa = -1.0 * Math.signum(_B.z()) / htf.R();
+ double kappa = -1.0 * Math.signum(B) / htf.R();
double theta = Math.PI / 2.0 - Math.atan(htf.slope());
perPar.setElement(0, 0, kappa);
perPar.setElement(0, 1, theta);
perPar.setElement(0, 2, htf.phi0());
- perPar.setElement(0, 3, htf.dca());
+ perPar.setElement(0, 3, -1.0 * htf.dca());
perPar.setElement(0, 4, htf.z0());
}
return perPar;
}
- public class PerigeeParams {
+ public static class PerigeeParams {
private final BasicMatrix _params;
- private PerigeeParams(HelicalTrackFit htf) {
- _params = GBLOutput.this.getPerParVector(htf);
+ public PerigeeParams(HelicalTrackFit htf, double B) {
+ _params = getPerParVector(htf, B);
}
public BasicMatrix getParams() {
@@ -867,7 +867,7 @@
* @param htf input helix to find the track direction
* @return 3x3 projection matrix
*/
- Hep3Matrix getPerToClPrj(HelicalTrackFit htf) {
+ 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));
@@ -890,17 +890,17 @@
}
- public class ClParams {
+ public static class ClParams {
private BasicMatrix _params = new BasicMatrix(1, 5);
- private ClParams(HelicalTrackFit htf) {
+ public ClParams(HelicalTrackFit htf, double B) {
if (htf == null) {
return;
}
- Hep3Matrix perToClPrj = GBLOutput.this.getPerToClPrj(htf);
+ Hep3Matrix perToClPrj = getPerToClPrj(htf);
double d0 = -1 * htf.dca(); //sign convention for curvilinear frame
double z0 = htf.z0();
@@ -915,7 +915,7 @@
double lambda = Math.atan(htf.slope());
double q = Math.signum(htf.R());
- double qOverP = q / htf.p(Math.abs(_B.z()));
+ double qOverP = q / htf.p(Math.abs(B));
double phi = htf.phi0();
_params.setElement(0, 0, qOverP);
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 Wed Sep 2 17:41:06 2015
@@ -1,12 +1,13 @@
package org.hps.recon.tracking.gbl;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import static java.lang.Math.abs;
import static java.lang.Math.abs;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.VecOp;
-
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
@@ -14,7 +15,6 @@
import java.util.logging.Formatter;
import java.util.logging.Level;
import java.util.logging.Logger;
-
import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.SymMatrix;
import org.hps.recon.tracking.gbl.matrix.Vector;
@@ -37,12 +37,13 @@
* @author Norman A Graf
* @author Per Hansson Adrian <[log in to unmask]>
*
- * @version $Id$
+ * @version $Id: HpsGblRefitter.java 3460 2015-08-29 01:45:39Z
+ * [log in to unmask] $
*/
-public class HpsGblRefitter extends Driver
-{
- static Formatter f = new BasicLogFormatter();
- private static Logger logger = LogUtil.create(HpsGblRefitter.class.getSimpleName(), f,Level.WARNING);
+public class HpsGblRefitter extends Driver {
+
+ static Formatter f = new BasicLogFormatter();
+ private static Logger logger = LogUtil.create(HpsGblRefitter.class.getSimpleName(), f, Level.WARNING);
//private static final Logger logger = Logger.getLogger(HpsGblRefitter.class.getName());
private boolean _debug = false;
private final String trackCollectionName = "MatchedTracks";
@@ -52,63 +53,55 @@
private MilleBinary mille;
private String milleBinaryFileName = MilleBinary.DEFAULT_OUTPUT_FILE_NAME;
private boolean writeMilleBinary = false;
-
+
private final MakeGblTracks _makeTracks;
- public void setDebug(boolean debug)
- {
+ public void setDebug(boolean debug) {
_debug = debug;
_makeTracks.setDebug(debug);
}
- public void setMilleBinaryFileName(String filename)
- {
+ public void setMilleBinaryFileName(String filename) {
milleBinaryFileName = filename;
}
- public void setWriteMilleBinary(boolean writeMillepedeFile)
- {
+ public void setWriteMilleBinary(boolean writeMillepedeFile) {
writeMilleBinary = writeMillepedeFile;
}
- public HpsGblRefitter()
- {
+ public HpsGblRefitter() {
_makeTracks = new MakeGblTracks();
_makeTracks.setDebug(_debug);
logger.setLevel(Level.WARNING);
System.out.println("level " + logger.getLevel().toString());
}
-
+
//@Override
//public void setLogLevel(String logLevel) {
// logger.setLevel(Level.parse(logLevel));
//}
-
@Override
- protected void startOfData()
- {
+ protected void startOfData() {
if (writeMilleBinary) {
mille = new MilleBinary(milleBinaryFileName);
}
}
@Override
- protected void endOfData()
- {
+ protected void endOfData() {
if (writeMilleBinary) {
mille.close();
}
}
@Override
- protected void process(EventHeader event)
- {
+ protected void process(EventHeader event) {
// Hep3Vector bfieldvec = event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0., 1.));
Hep3Vector bfieldvec = event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0., 500.));//mg...get the b-field in the middle of the magnet
double bfield = bfieldvec.y();
//double bfac = 0.0002998 * bfield;
double bfac = Constants.fieldConversion * bfield;
-
+
// get the tracks
if (!event.hasCollection(Track.class, trackCollectionName)) {
if (_debug) {
@@ -128,13 +121,12 @@
return;
}
-
List<LCRelation> track2GblTrackRelations = event.get(LCRelation.class, track2GblTrackRelationName);
//need a map of GBLTrackData keyed on the Generic object from which it created
Map<GenericObject, GBLTrackData> gblObjMap = new HashMap<GenericObject, GBLTrackData>();
//need a map of SeedTrack to GBLTrackData keyed on the track object from which it created
- Map<GBLTrackData, Track> gblToSeedMap = new HashMap<GBLTrackData, Track>();
-
+ Map<GBLTrackData, Track> gblToSeedMap = new HashMap<GBLTrackData, Track>();
+
// loop over the relations
for (LCRelation relation : track2GblTrackRelations) {
Track t = (Track) relation.getFrom();
@@ -162,57 +154,51 @@
stripsGblMap.get(gblT).add(sd);
}
}
-
- // loop over the tracks and do the GBL fit
+
+ // loop over the tracks and do the GBL fit
List<FittedGblTrajectory> trackFits = new ArrayList<FittedGblTrajectory>();
int trackNum = 0;
logger.info("Trying to fit " + stripsGblMap.size() + " tracks");
for (GBLTrackData t : stripsGblMap.keySet()) {
FittedGblTrajectory traj = fit(stripsGblMap.get(t), bfac);
++trackNum;
- if(traj!=null) {
+ if (traj != null) {
logger.info("GBL fit successful");
- if(_debug) System.out.printf("%s: GBL fit successful.\n",getClass().getSimpleName());
+ if (_debug) {
+ System.out.printf("%s: GBL fit successful.\n", getClass().getSimpleName());
+ }
traj.set_seed(gblToSeedMap.get(t));
traj.set_track_data(t);
trackFits.add(traj);
} else {
logger.info("GBL fit failed");
- if(_debug) System.out.printf("%s: GBL fit failed.\n",getClass().getSimpleName());
- }
- }
-
+ if (_debug) {
+ System.out.printf("%s: GBL fit failed.\n", getClass().getSimpleName());
+ }
+ }
+ }
+
logger.info(event.get(Track.class, trackCollectionName).size() + " tracks in collection \"" + trackCollectionName + "\"");
logger.info(gblObjMap.size() + " tracks in gblObjMap");
logger.info(gblToSeedMap.size() + " tracks in gblToSeedMap");
logger.info(stripsGblMap.size() + " tracks in stripsGblMap");
logger.info(trackFits.size() + " fitted GBL tracks before adding to event");
-
-
-
-
+
_makeTracks.Process(event, trackFits, bfield);
- if(_debug) System.out.printf("%s: Done.\n",getClass().getSimpleName());
-
-
- }
-
- private FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac)
- {
+ if (_debug) {
+ System.out.printf("%s: Done.\n", getClass().getSimpleName());
+ }
+
+ }
+
+ private FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac) {
// path length along trajectory
double s = 0.;
// jacobian to transport errors between points along the path
Matrix jacPointToPoint = new Matrix(5, 5);
jacPointToPoint.UnitMatrix();
- // Option to use uncorrelated MS errors
- // This is similar to what is done in lcsim seedtracker
- // The msCov below holds the MS errors
- // This is for testing purposes only.
- boolean useUncorrMS = false;
- Matrix msCov = new Matrix(5, 5);
- Matrix measMsCov = new Matrix(2, 2);
// Vector of the strip clusters used for the GBL fit
List<GblPoint> listOfPoints = new ArrayList<GblPoint>();
@@ -236,13 +222,12 @@
if (_debug) {
System.out.println("HpsGblFitter: " + "Path length step " + step + " from " + s + " to " + strip.getPath3D());
}
-
+
// get measurement frame unit vectors
Hep3Vector u = strip.getU();
Hep3Vector v = strip.getV();
Hep3Vector w = strip.getW();
-
-
+
// Measurement direction (perpendicular and parallel to strip direction)
Matrix mDir = new Matrix(2, 3);
mDir.set(0, 0, u.x());
@@ -327,7 +312,7 @@
System.out.println("HpsGblFitter: " + "measPrec:");
measPrec.print(4, 6);
}
-
+
//Find the Jacobian to be able to propagate the covariance matrix to this strip position
jacPointToPoint = gblSimpleJacobianLambdaPhi(step, cosLambda, abs(bfac));
@@ -338,6 +323,13 @@
// Get the transpose of the Jacobian
Matrix jacPointToPointTransposed = jacPointToPoint.copy().transpose();
+
+ // Option to use uncorrelated MS errors
+ // This is similar to what is done in lcsim seedtracker
+ // The msCov below holds the MS errors
+ // This is for testing purposes only.
+ boolean useUncorrMS = false;
+ Matrix msCov = new Matrix(5, 5);
// Propagate the MS covariance matrix (in the curvilinear frame) to this strip position
msCov = msCov.times(jacPointToPointTransposed);
@@ -346,7 +338,7 @@
// Get the MS covariance for the measurements in the measurement frame
Matrix proL2mTransposed = proL2m.copy().transpose();
- measMsCov = proL2m.times((msCov.getMatrix(3, 4, 3, 4)).times(proL2mTransposed));
+ Matrix measMsCov = proL2m.times((msCov.getMatrix(3, 4, 3, 4)).times(proL2mTransposed));
if (_debug) {
System.out.println("HpsGblFitter: " + " msCov at this point:");
@@ -388,54 +380,51 @@
msCov.set(1, 1, msCov.get(1, 1) + scatErr.get(0) * scatErr.get(0));
msCov.set(2, 2, msCov.get(2, 2) + scatErr.get(1) * scatErr.get(1));
-
-
- //// Calculate global derivatives for this point
-
+ //// Calculate global derivatives for this point
// track direction in tracking/global frame
- Hep3Vector tDirGlobal = new BasicHep3Vector(cosPhi * cosLambda, sinPhi * cosLambda, sinLambda);
-
+ Hep3Vector tDirGlobal = new BasicHep3Vector(cosPhi * cosLambda, sinPhi * cosLambda, sinLambda);
+
// Cross-check that the input is consistent
- if( VecOp.sub(tDirGlobal, strip.getTrackDirection()).magnitude() > 0.00001) {
- throw new RuntimeException("track directions are inconsistent: " + tDirGlobal.toString() + " and " + strip.getTrackDirection().toString());
- }
- // rotate track direction to measurement frame
- Hep3Vector tDirMeas = new BasicHep3Vector(VecOp.dot(tDirGlobal, u), VecOp.dot(tDirGlobal, v), VecOp.dot(tDirGlobal, w));
- // TODO this is a trivial one. Fix it.
- Hep3Vector normalMeas = new BasicHep3Vector(VecOp.dot(w, u), VecOp.dot(w, v), VecOp.dot(w, w));
-
- // vector coplanar with measurement plane from origin to prediction
- Hep3Vector tPosMeas = strip.getTrackPos();
-
- // measurements: non-measured directions
- double vmeas = 0.;
- double wmeas = 0.;
-
- // calculate and add derivatives to point
- GlobalDers glDers = new GlobalDers(strip.getId(), meas.get(0), vmeas, wmeas, tDirMeas, tPosMeas, normalMeas);
-
- //TODO find a more robust way to get half.
- boolean isTop = Math.sin(strip.getTrackLambda()) > 0 ? true : false;
-
- // Get the list of millepede parameters
- List<MilleParameter> milleParameters = glDers.getDers(isTop);
-
- // need to make vector and matrices for interface
- List<Integer> labGlobal = new ArrayList<Integer>();
- Matrix addDer = new Matrix(1, milleParameters.size());
- for(int i=0; i < milleParameters.size(); ++i) {
- labGlobal.add(milleParameters.get(i).getId());
- addDer.set(0, i, milleParameters.get(i).getValue());
- }
- point.addGlobals(labGlobal, addDer);
- String logders = "";
- for(int i=0; i < milleParameters.size(); ++i) {
- logders += labGlobal.get(i) + "\t" + addDer.get(0, i) + "\n";
- }
- logger.info("\n"+ logders);
-
- logger.info("uRes " + strip.getId() + " uRes " + uRes + " pred (" + strip.getTrackPos().x() + "," + strip.getTrackPos().y() + "," + strip.getTrackPos().z() + ") s(3D) " + strip.getPath3D());
-
+ if (VecOp.sub(tDirGlobal, strip.getTrackDirection()).magnitude() > 0.00001) {
+ throw new RuntimeException("track directions are inconsistent: " + tDirGlobal.toString() + " and " + strip.getTrackDirection().toString());
+ }
+ // rotate track direction to measurement frame
+ Hep3Vector tDirMeas = new BasicHep3Vector(VecOp.dot(tDirGlobal, u), VecOp.dot(tDirGlobal, v), VecOp.dot(tDirGlobal, w));
+ // TODO this is a trivial one. Fix it.
+ Hep3Vector normalMeas = new BasicHep3Vector(VecOp.dot(w, u), VecOp.dot(w, v), VecOp.dot(w, w));
+
+ // vector coplanar with measurement plane from origin to prediction
+ Hep3Vector tPosMeas = strip.getTrackPos();
+
+ // measurements: non-measured directions
+ double vmeas = 0.;
+ double wmeas = 0.;
+
+ // calculate and add derivatives to point
+ GlobalDers glDers = new GlobalDers(strip.getId(), meas.get(0), vmeas, wmeas, tDirMeas, tPosMeas, normalMeas);
+
+ //TODO find a more robust way to get half.
+ boolean isTop = Math.sin(strip.getTrackLambda()) > 0;
+
+ // Get the list of millepede parameters
+ List<MilleParameter> milleParameters = glDers.getDers(isTop);
+
+ // need to make vector and matrices for interface
+ List<Integer> labGlobal = new ArrayList<Integer>();
+ Matrix addDer = new Matrix(1, milleParameters.size());
+ for (int i = 0; i < milleParameters.size(); ++i) {
+ labGlobal.add(milleParameters.get(i).getId());
+ addDer.set(0, i, milleParameters.get(i).getValue());
+ }
+ point.addGlobals(labGlobal, addDer);
+ String logders = "";
+ for (int i = 0; i < milleParameters.size(); ++i) {
+ logders += labGlobal.get(i) + "\t" + addDer.get(0, i) + "\n";
+ }
+ logger.info("\n" + logders);
+
+ logger.info("uRes " + strip.getId() + " uRes " + uRes + " pred (" + strip.getTrackPos().x() + "," + strip.getTrackPos().y() + "," + strip.getTrackPos().z() + ") s(3D) " + strip.getPath3D());
+
//go to next point
s += step;
@@ -460,7 +449,7 @@
double[] dVals = new double[2];
int[] iVals = new int[1];
traj.fit(dVals, iVals, "");
- logger.info("fit result: Chi2="+ dVals[0] + " Ndf=" + iVals[0] + " Lost=" + dVals[1]);
+ logger.info("fit result: Chi2=" + dVals[0] + " Ndf=" + iVals[0] + " Lost=" + dVals[1]);
Vector aCorrection = new Vector(5);
SymMatrix aCovariance = new SymMatrix(5);
traj.getResults(1, aCorrection, aCovariance);
@@ -472,7 +461,7 @@
}
logger.fine("locPar " + aCorrection.toString());
-
+
// // write to MP binary file
if (writeMilleBinary) {
traj.milleOut(mille);
@@ -482,12 +471,10 @@
}
@Override
- protected void detectorChanged(Detector detector)
- {
- }
-
- private Matrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac)
- {
+ protected void detectorChanged(Detector detector) {
+ }
+
+ private Matrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac) {
/**
* Simple jacobian: quadratic in arc length difference. using lambda phi
* as directions
@@ -514,21 +501,21 @@
jac.set(4, 1, ds);
return jac;
}
-
-
+
private static class GlobalDers {
- private int _layer;
- private double _umeas; // measurement direction
- private double _vmeas; // unmeasured direction
- private double _wmeas; // normal to plane
- private Hep3Vector _t; // track direction
- private Hep3Vector _p; // track prediction
- private Hep3Vector _n; // normal to plane
- private Matrix _dm_dg; //Global derivaties of the local measurements
- private Matrix _dr_dm; //Derivatives of residuals w.r.t. measurement
- private Matrix _dr_dg; //Derivatives of residuals w.r.t. global parameters
-
- public GlobalDers(int layer,double umeas, double vmeas, double wmeas, Hep3Vector tDir, Hep3Vector tPred, Hep3Vector normal) {
+
+ private final int _layer;
+ private final double _umeas; // measurement direction
+ private final double _vmeas; // unmeasured direction
+ private final double _wmeas; // normal to plane
+ private final Hep3Vector _t; // track direction
+ private final Hep3Vector _p; // track prediction
+ private final Hep3Vector _n; // normal to plane
+ private final Matrix _dm_dg; //Global derivaties of the local measurements
+ private final Matrix _dr_dm; //Derivatives of residuals w.r.t. measurement
+ private final Matrix _dr_dg; //Derivatives of residuals w.r.t. global parameters
+
+ public GlobalDers(int layer, double umeas, double vmeas, double wmeas, Hep3Vector tDir, Hep3Vector tPred, Hep3Vector normal) {
_layer = layer;
_umeas = umeas;
_vmeas = vmeas;
@@ -541,7 +528,7 @@
// Derivatives of perturbed measurement w.r.t. global parameters
_dm_dg = getMeasDers();
// Calculate, by chain rule, derivatives of residuals w.r.t. global parameters
- _dr_dg = _dr_dm.times(_dm_dg);
+ _dr_dg = _dr_dm.times(_dm_dg);
//logger.log(Level.FINER," dr_dm\n"+ _dr_dm.toString() + "\ndm_dg\n" + _dm_dg.toString() + "\ndr_dg\n" +_dr_dg.toString());
//logger.info("loglevel " + logger.getLevel().toString());
@@ -553,9 +540,9 @@
//print self.dr_dg
}
-
/**
- * Derivative of mt, the perturbed measured coordinate vector m w.r.t. to global parameters: u,v,w,alpha,beta,gamma
+ * Derivative of mt, the perturbed measured coordinate vector m w.r.t.
+ * to global parameters: u,v,w,alpha,beta,gamma
*/
private Matrix getMeasDers() {
@@ -608,27 +595,28 @@
}
/**
- * Derivatives of the local perturbed residual w.r.t. the measurements m (u,v,w)'
+ * Derivatives of the local perturbed residual w.r.t. the measurements m
+ * (u,v,w)'
*/
private Matrix getResDers() {
- double tdotn = VecOp.dot(_t, _n);
- Matrix dr_dm = Matrix.identity(3,3);
+ double tdotn = VecOp.dot(_t, _n);
+ Matrix dr_dm = Matrix.identity(3, 3);
//print 't ', self.t, ' n ', self.n, ' dot(t,n) ', tdotn
//logger.info("t " + _t.toString() +" n " + _n.toString() + " dot(t,n) " + tdotn);
double delta, val;
- for(int i=0; i<3; ++i) {
- for(int j=0; j<3; ++j) {
- delta = i==j ? 1. : 0.;
- val = delta - _t.v()[i] * _n.v()[j]/tdotn;
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ delta = i == j ? 1. : 0.;
+ val = delta - _t.v()[i] * _n.v()[j] / tdotn;
dr_dm.set(i, j, val);
}
}
return dr_dm;
}
-
/**
* Turn derivative matrix into @Milleparameter
+ *
* @param isTop - top or bottom track
* @return list of @Milleparameters
*/
@@ -639,16 +627,16 @@
double value;
List<MilleParameter> milleParameters = new ArrayList<MilleParameter>();
int topBot = isTop == true ? 1 : 2;
- for(int ip=1; ip<7; ++ip) {
- if(ip > 3) {
+ for (int ip = 1; ip < 7; ++ip) {
+ if (ip > 3) {
transRot = 2;
- direction = ((ip-1) % 3) + 1;
+ direction = ((ip - 1) % 3) + 1;
} else {
transRot = 1;
direction = ip;
}
label = topBot * MilleParameter.half_offset + transRot * MilleParameter.type_offset + direction * MilleParameter.dimension_offset + _layer;
- value = _dr_dg.get(0, ip-1);
+ value = _dr_dg.get(0, ip - 1);
MilleParameter milleParameter = new MilleParameter(label, value, 0.0);
milleParameters.add(milleParameter);
}
@@ -658,114 +646,111 @@
/*
- class globalDers:
- def __init__(self,layer,umeas,vmeas,wmeas, tDir, tPred, normal):
- self.layer = layer # measurement direction
- self.umeas = umeas # measurement direction
- self.vmeas = vmeas # unmeasured direction
- self.wmeas = wmeas # normal to plane
- self.t = tDir # track direction
- self.p = tPred # track prediction
- self.n = normal # normal to plane
- # Global derivaties of the local measurements
- self.dm_dg = self.getMeasDers()
- # Derivatives of residuals w.r.t. measurement
- self.dr_dm = self.getResDers()
- # Derivatives of residuals w.r.t. global parameters
- self.dr_dg = np.dot(self.dr_dm, self.dm_dg)
- #print 'dm_dg'
- #print dm_dg
- #print 'dr_dm'
- #print dr_dm
- #print 'dr_dg'
- #print self.dr_dg
-
- def dump(self):
- print 'globalDers:'
- print 'layer ', self.layer
- print 'umeas ', self.umeas, ' vmeas ', self.vmeas, ' wmeas ', self.wmeas
- print 't ', self.t, ' p ', self.p, ' n ', self.n
- print 'dm_dg\n',self.dm_dg, '\ndr_dm\n',self.dr_dm,'\ndr_dg\n',self.dr_dg
-
- def getDers(self,isTop):
- half_offset = 10000
- translation_offset = 1000
- direction_offset = 100
- topBot = 1
- transRot = 1
- direction = 1
- if not isTop:
- topBot = 2
- res = {}
- labels = []
- ders = []
- for ip, name in global_params.iteritems():
- if ip > 3:
- transRot = 2
- direction = ((ip-1) % 3) + 1
- else:
- direction = ip
- label = (int)(topBot * half_offset + transRot * translation_offset + direction * direction_offset + self.layer)
- labels.append(label)
- ders.append(self.dr_dg[0,ip-1])
-
- return {'labels':np.array([labels]) , 'ders':np.array([ders])}
-
-
-
-
- def getResDers(self):
- # Derivatives of the local perturbed residual w.r.t. the measurements m (u,v,w)'
- tdotn = np.dot(self.t.T,self.n)
- drdg = np.eye(3)
- #print 't ', self.t, ' n ', self.n, ' dot(t,n) ', tdotn
- for i in range(3):
- for j in range(3):
- delta = 0.
- if i==j:
- delta = 1.
- drdg[i][j] = delta - self.t[i]*self.n[j]/tdotn[0]
- return drdg
-
-
-
-
- def getMeasDers(self):
- # Derivative of mt, the perturbed measured coordinate vector m
- # w.r.t. to global parameters: u,v,w,alpha,beta,gamma
-
- # Derivative of the local measurement for a translation in u
- dmu_du = 1.
- dmv_du = 0.
- dmw_du = 0.
- # Derivative of the local measurement for a translation in v
- dmu_dv = 0.
- dmv_dv = 1.
- dmw_dv = 0.
- # Derivative of the local measurement for a translation in w
- dmu_dw = 0.
- dmv_dw = 0.
- dmw_dw = 1.
- # Derivative of the local measurement for a rotation around u-axis (alpha)
- dmu_dalpha = 0.
- dmv_dalpha = self.p[2] # self.wmeas
- dmw_dalpha = -1.0 * self.p[1] # -1.0 * self.vmeas
- # Derivative of the local measurement for a rotation around v-axis (beta)
- dmu_dbeta = -1.0 * self.p[2] #-1.0 * self.wmeas
- dmv_dbeta = 0.
- dmw_dbeta = self.p[0] #self.umeas
- # Derivative of the local measurement for a rotation around w-axis (gamma)
- dmu_dgamma = self.p[1] # self.vmeas
- dmv_dgamma = -1.0 * self.p[0] # -1.0 * self.umeas
- dmw_dgamma = 0.
- # put into matrix
- dmdg = np.array([[dmu_du, dmu_dv, dmu_dw, dmu_dalpha, dmu_dbeta, dmu_dgamma],[dmv_du, dmv_dv, dmv_dw, dmv_dalpha, dmv_dbeta, dmv_dgamma],[dmw_du, dmw_dv, dmw_dw, dmw_dalpha, dmw_dbeta, dmw_dgamma]])
- #print dmw_dbeta
- #dmdg = np.array([[dmu_du, dmu_dv],[dmu_dw, dmu_dalpha], [dmw_dbeta, dmw_dgamma]])
- return dmdg
+ class globalDers:
+ def __init__(self,layer,umeas,vmeas,wmeas, tDir, tPred, normal):
+ self.layer = layer # measurement direction
+ self.umeas = umeas # measurement direction
+ self.vmeas = vmeas # unmeasured direction
+ self.wmeas = wmeas # normal to plane
+ self.t = tDir # track direction
+ self.p = tPred # track prediction
+ self.n = normal # normal to plane
+ # Global derivaties of the local measurements
+ self.dm_dg = self.getMeasDers()
+ # Derivatives of residuals w.r.t. measurement
+ self.dr_dm = self.getResDers()
+ # Derivatives of residuals w.r.t. global parameters
+ self.dr_dg = np.dot(self.dr_dm, self.dm_dg)
+ #print 'dm_dg'
+ #print dm_dg
+ #print 'dr_dm'
+ #print dr_dm
+ #print 'dr_dg'
+ #print self.dr_dg
+
+ def dump(self):
+ print 'globalDers:'
+ print 'layer ', self.layer
+ print 'umeas ', self.umeas, ' vmeas ', self.vmeas, ' wmeas ', self.wmeas
+ print 't ', self.t, ' p ', self.p, ' n ', self.n
+ print 'dm_dg\n',self.dm_dg, '\ndr_dm\n',self.dr_dm,'\ndr_dg\n',self.dr_dg
+
+ def getDers(self,isTop):
+ half_offset = 10000
+ translation_offset = 1000
+ direction_offset = 100
+ topBot = 1
+ transRot = 1
+ direction = 1
+ if not isTop:
+ topBot = 2
+ res = {}
+ labels = []
+ ders = []
+ for ip, name in global_params.iteritems():
+ if ip > 3:
+ transRot = 2
+ direction = ((ip-1) % 3) + 1
+ else:
+ direction = ip
+ label = (int)(topBot * half_offset + transRot * translation_offset + direction * direction_offset + self.layer)
+ labels.append(label)
+ ders.append(self.dr_dg[0,ip-1])
+
+ return {'labels':np.array([labels]) , 'ders':np.array([ders])}
+
+
+
+
+ def getResDers(self):
+ # Derivatives of the local perturbed residual w.r.t. the measurements m (u,v,w)'
+ tdotn = np.dot(self.t.T,self.n)
+ drdg = np.eye(3)
+ #print 't ', self.t, ' n ', self.n, ' dot(t,n) ', tdotn
+ for i in range(3):
+ for j in range(3):
+ delta = 0.
+ if i==j:
+ delta = 1.
+ drdg[i][j] = delta - self.t[i]*self.n[j]/tdotn[0]
+ return drdg
+
+
+
+
+ def getMeasDers(self):
+ # Derivative of mt, the perturbed measured coordinate vector m
+ # w.r.t. to global parameters: u,v,w,alpha,beta,gamma
+
+ # Derivative of the local measurement for a translation in u
+ dmu_du = 1.
+ dmv_du = 0.
+ dmw_du = 0.
+ # Derivative of the local measurement for a translation in v
+ dmu_dv = 0.
+ dmv_dv = 1.
+ dmw_dv = 0.
+ # Derivative of the local measurement for a translation in w
+ dmu_dw = 0.
+ dmv_dw = 0.
+ dmw_dw = 1.
+ # Derivative of the local measurement for a rotation around u-axis (alpha)
+ dmu_dalpha = 0.
+ dmv_dalpha = self.p[2] # self.wmeas
+ dmw_dalpha = -1.0 * self.p[1] # -1.0 * self.vmeas
+ # Derivative of the local measurement for a rotation around v-axis (beta)
+ dmu_dbeta = -1.0 * self.p[2] #-1.0 * self.wmeas
+ dmv_dbeta = 0.
+ dmw_dbeta = self.p[0] #self.umeas
+ # Derivative of the local measurement for a rotation around w-axis (gamma)
+ dmu_dgamma = self.p[1] # self.vmeas
+ dmv_dgamma = -1.0 * self.p[0] # -1.0 * self.umeas
+ dmw_dgamma = 0.
+ # put into matrix
+ dmdg = np.array([[dmu_du, dmu_dv, dmu_dw, dmu_dalpha, dmu_dbeta, dmu_dgamma],[dmv_du, dmv_dv, dmv_dw, dmv_dalpha, dmv_dbeta, dmv_dgamma],[dmw_du, dmw_dv, dmw_dw, dmw_dalpha, dmw_dbeta, dmw_dgamma]])
+ #print dmw_dbeta
+ #dmdg = np.array([[dmu_du, dmu_dv],[dmu_dw, dmu_dalpha], [dmw_dbeta, dmw_dgamma]])
+ return dmdg
*/
}
}
-
-
-
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 Wed Sep 2 17:41:06 2015
@@ -5,12 +5,14 @@
import hep.physics.vec.Hep3Matrix;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
-
import java.util.ArrayList;
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.gbl.GBLOutput.ClParams;
+import org.hps.recon.tracking.gbl.GBLOutput.PerigeeParams;
+import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.SymMatrix;
import org.hps.recon.tracking.gbl.matrix.Vector;
import org.hps.util.BasicLogFormatter;
@@ -93,11 +95,9 @@
// Retrieve the helix and save the relevant bits of helix info
HelicalTrackFit helix = trackseed.getHelix();
- double gblParameters[] = getGblCorrectedHelixParameters(helix, fittedTraj, bfield);
- trk.setTrackParameters(gblParameters, bfield); // Sets first TrackState.
- //TODO Use GBL covariance matrix
- //SymmetricMatrix gblCovariance = getGblCorrectedHelixCovariance(helix, fittedTraj, bfield);
- trk.setCovarianceMatrix(helix.covariance()); // Modifies first TrackState.
+ Pair<double[], SymmetricMatrix> correctedHelixParams = getGblCorrectedHelixParameters(helix, fittedTraj, bfield);
+ trk.setTrackParameters(correctedHelixParams.getFirst(), bfield); // Sets first TrackState.
+ trk.setCovarianceMatrix(correctedHelixParams.getSecond()); // Modifies first TrackState.
trk.setChisq(fittedTraj.get_chi2());
trk.setNDF(fittedTraj.get_ndf());
@@ -138,49 +138,44 @@
}
/**
- * Compute the track fit covariance matrix
- *
- * @param helix - original seed track
- * @param traj - fitted GBL trajectory
- * @return covariance matrix.
- */
- private SymmetricMatrix getGblCorrectedHelixCovariance(
- HelicalTrackFit helix, FittedGblTrajectory traj, double bfield) {
- // TODO Actually implement this method
- return helix.covariance();
- }
-
- /**
- * Compute the updated helix parameters.
+ * Compute the updated helix parameters and covariance matrix.
*
* @param helix - original seed track
* @param traj - fitted GBL trajectory
* @return corrected parameters.
*/
- private double[] getGblCorrectedHelixParameters(HelicalTrackFit helix, FittedGblTrajectory traj, double bfield) {
+ private Pair<double[], SymmetricMatrix> getGblCorrectedHelixParameters(HelicalTrackFit helix, FittedGblTrajectory traj, double bfield) {
// get seed helix parameters
- double d0 = -1.0 * helix.dca(); // correct for different sign convention of d0 in curvilinear frame
+// double p = helix.p(Math.abs(bfield));
+// double q = Math.signum(helix.R());
+// double qOverP = q / p;
+// ClParams clParams = new ClParams(helix, bfield);
+// PerigeeParams perParams = new PerigeeParams(helix, bfield);
+ 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());
- double p = helix.p(Math.abs(bfield));
- double q = Math.signum(helix.R());
- double qOverP = q / p;
-
+
+// System.out.println("clParams: " + clParams.getParams());
+// System.out.println("perParams: " + perParams.getParams());
+// System.out.format("converted params: qOverP %f, d0 %f, z0 %f, phi0 %f, lambda %f\n", qOverP, d0, z0, phi0, lambda);
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
Vector locPar = new Vector(5);
SymMatrix locCov = new SymMatrix(5);
traj.get_traj().getResults(1, 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());
- double xTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
- double yTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
-
- logger.info((helix.slope() > 0 ? "top: " : "bot ") + "qOverPCorr " + qOverPCorr + " xTCorr " + xTCorr + " yTCorr " + yTCorr + " xtPrimeCorr " + xTPrimeCorr + " yTPrimeCorr " + yTPrimeCorr);
+
+ 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();
@@ -204,22 +199,43 @@
// calculate new curvature
double qOverP_gbl = qOverP + qOverPCorr;
- double pt_gbl = Math.abs(1.0 / qOverP_gbl) * Math.sin((Math.PI / 2.0 - lambda_gbl));
- double C_gbl = Constants.fieldConversion * bfield / pt_gbl;
- //make sure sign is not changed
- C_gbl = Math.signum(qOverP_gbl) * Math.abs(C_gbl);
+// 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);
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)));
+
+ 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.z0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), clToPerPrj.e(2, 0));
+ jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(2, 1));
+ jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1.0);
+ jacobian.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Math.pow(Math.cos(lambda_gbl), -2.0));
+ 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));
+
+ 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.curvatureIndex] = C_gbl;
parameters_gbl[HelicalTrackFit.slopeIndex] = slope_gbl;
- parameters_gbl[HelicalTrackFit.phi0Index] = phi0_gbl;
- return parameters_gbl;
+
+ return new Pair<double[], SymmetricMatrix>(parameters_gbl, cov);
}
public void setTrkCollectionName(String name) {
|
