Author: [log in to unmask]
Date: Thu Sep 3 09:55:45 2015
New Revision: 3505
Log:
calculate covariance matrix template
Modified:
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/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 Thu Sep 3 09:55:45 2015
@@ -1,6 +1,7 @@
package org.hps.recon.tracking.gbl;
import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Matrix;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Matrix;
import hep.physics.vec.Hep3Vector;
@@ -26,6 +27,7 @@
import org.lcsim.event.TrackerHit;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelixUtils;
import org.lcsim.lcio.LCIOConstants;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -216,6 +218,50 @@
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 SymMatrix(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,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));
+ 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.curvatureIndex,FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), -1*Bz/cosLambda);
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex,FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1*q*Bz*tanLambda/(p*cosLambda));
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex,FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 0);
+ 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.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.set(HelicalTrackFit.phi0Index,FittedGblTrajectory.GBLPARIDX.YTPRIME.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.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), VdotK/TdotI);
+ covariance_gbl.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -1/TdotI);
+
+
+
+ // 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));
