lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.10 -r1.11
--- HitUtils.java 21 May 2009 21:48:15 -0000 1.10
+++ HitUtils.java 16 Nov 2010 16:23:18 -0000 1.11
@@ -173,21 +173,28 @@
double separation = SensorSeperation(strip1, strip2);
// Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
double salpha = v1Dotu2(strip1, strip2);
- // Calculate the scale factor: factor = (1 + gamma) / 2 * sin(alpha)
- double factor = (1 + gamma) / (2. * salpha);
- // Calculate the covariance matrix from resolution: cov = factor * (v2 * v2^T * du1^2 + v1 * v1^T * du2^2)
- BasicMatrix cov = Vec2Matrix(VecOp.mult(strip1.du() * factor, strip2.v()));
- MatrixAdd(Vec2Matrix(VecOp.mult(strip2.du() * factor, strip1.v())), cov);
+ // Calculate the scale factor: factor = (1 + gamma)^2 / 4 * sin^2(alpha)
+ double factor = (1 + gamma)*(1 + gamma) / (4. * salpha*salpha);
+ // Calculate v * v^T for strips 1 and 2
+ Matrix v1 = v2m(strip1.v());
+ Matrix v2 = v2m(strip2.v());
+ Matrix v1v1t = MatrixOp.mult(v1, MatrixOp.transposed(v1));
+ Matrix v2v2t = MatrixOp.mult(v2, MatrixOp.transposed(v2));
+ // Find measurement uncertainties for strips 1 and 2
+ double du1 = strip1.du();
+ double du2 = strip2.du();
// Calculate the uncertainty in the unmeasured coordinate due to not knowing the track direction
// by assuming phat . u has an uncertainty of 2/sqrt(12) so dv = 2 / sqrt(12) * seperation / sin(alpha)
double dv = Math.abs(2. * separation / (Math.sqrt(12) * salpha));
// Don't let dv by greater than the strip length / sqrt(12)
double dv1 = Math.min(dv, (strip1.vmax()-strip1.vmin()) / Math.sqrt(12.));
double dv2 = Math.min(dv, (strip2.vmax()-strip2.vmin()) / Math.sqrt(12.));
- // Add the covariance contributions from not knowing the track direction
- // cov += v1 * v1^T * (dv1/2)^2 + v2 * v2^T * (dv2/2)^2
- MatrixAdd(Vec2Matrix(VecOp.mult(0.5 * dv1, strip1.v())), cov);
- MatrixAdd(Vec2Matrix(VecOp.mult(0.5 * dv2, strip2.v())), cov);
+ // Calculate the covariance matrix.
+ // From resolution: cov = factor * (v2 * v2^T * du1^2 + v1 * v1^T * du2^2)
+ // From direction: + (v1 * v1^T * (dv1/2)^2 + v2 * v2^T * (dv2/2)^2
+ Matrix cov1 = MatrixOp.mult(factor * du2*du2 + 0.25 * dv1*dv1, v1v1t);
+ Matrix cov2 = MatrixOp.mult(factor * du1*du1 + 0.25 * dv2*dv2, v2v2t);
+ Matrix cov = MatrixOp.add(cov1, cov2);
return new SymmetricMatrix(cov);
}
@@ -229,9 +236,9 @@
Hep3Vector dir = trkdir.Direction();
Matrix dirderiv = trkdir.Derivatives();
// Calculate phat x v1
- Hep3Vector pcrossv1 = VecOp.cross(dir, strip1.v());
+ Matrix pcrossv1 = v2m(VecOp.cross(dir, strip1.v()));
// Calculate phat x v2
- Hep3Vector pcrossv2 = VecOp.cross(dir, strip2.v());
+ Matrix pcrossv2 = v2m(VecOp.cross(dir, strip2.v()));
// Calculate phat . w
double pdotw = NonZeroDotProduct(dir, strip1.w());
// salpha is the sin of the stereo angle
@@ -240,16 +247,22 @@
double factor = SensorSeperation(strip1, strip2) / (2 * salpha * pdotw*pdotw);
// Create the matrix of position derivatives: d_i,j = dr_i / dphat_j
// The matrix d is given by factor * (v1 * (phat x v2)^T + v2 * (phat x v1)^T where ^T means transpose
- MutableMatrix d = Vec2Matrix(strip1.v(), VecOp.mult(factor, pcrossv2));
- MatrixAdd(Vec2Matrix(strip2.v(), VecOp.mult(factor, pcrossv1)), d);
+ Matrix v1 = v2m(strip1.v());
+ Matrix v2 = v2m(strip2.v());
+ Matrix d = MatrixOp.mult(factor, MatrixOp.add(MatrixOp.mult(v1, MatrixOp.transposed(pcrossv2)),
+ MatrixOp.mult(v2, MatrixOp.transposed(pcrossv1))));
Matrix dh = MatrixOp.mult(d, dirderiv);
// Construct the transpose of dh
- Matrix dht = MatrixTranspose(dh);
+ Matrix dht = MatrixOp.transposed(dh);
// Calculate the covariance contributions from the direction uncertainty: cov = dh * hcov * dh^T
- MutableMatrix cov = (MutableMatrix) MatrixOp.mult(dh, MatrixOp.mult(hcov, dht));
- // Add the contributions from measurement errors: cov += (v2 * v2^T * du1^2 + v1 * v1^T * du2^2) / sin(alpha)^2
- MatrixAdd(Vec2Matrix(VecOp.mult(strip1.du()/salpha, strip2.v())), cov);
- MatrixAdd(Vec2Matrix(VecOp.mult(strip2.du()/salpha, strip1.v())), cov);
+ Matrix cov_dir = MatrixOp.mult(dh, MatrixOp.mult(hcov, dht));
+ // Calculate the contributions from measurement errors: cov += (v2 * v2^T * du1^2 + v1 * v1^T * du2^2) / sin(alpha)^2
+ double du1 = strip1.du();
+ double du2 = strip2.du();
+ Matrix cov1 = MatrixOp.mult(du1*du1 / (salpha*salpha), MatrixOp.mult(v2, MatrixOp.transposed(v2)));
+ Matrix cov2 = MatrixOp.mult(du2*du2 / (salpha*salpha), MatrixOp.mult(v1, MatrixOp.transposed(v1)));
+ // Sum all the contributions
+ Matrix cov = MatrixOp.add(cov_dir, MatrixOp.add(cov1, cov2));
// Convert to a symmetric matrix
return new SymmetricMatrix(cov);
@@ -298,7 +311,7 @@
// Construct the matrix dh = d^T * dirderiv
Matrix dh = MatrixOp.mult(dT, dirderiv);
// Construct the transpose of dh
- Matrix dhT = MatrixTranspose(dh);
+ Matrix dhT = MatrixOp.transposed(dh);
// Calculate the uncertainty in v1 from the direction uncertainty: cov = dh * hcov * dh^T
MutableMatrix cov = (MutableMatrix) MatrixOp.mult(dh, MatrixOp.mult(hcov, dhT));
// Return the uncertainty in v1: dv1^2 = ((u1 . u2)^2 * du1^2 + du2^2) / sin^2(alpha) + cov
@@ -324,66 +337,12 @@
}
return cth;
}
-
- /**
- * Form the outer matrix product vec1 * vec2^T from the vectors
- * vec1 and vec2.
- * @param vec1 first vector
- * @param vec2 second vector
- * @return outer product
- */
- private static BasicMatrix Vec2Matrix(Hep3Vector vec1, Hep3Vector vec2) {
- BasicMatrix mat = new BasicMatrix(3, 3);
- for (int i = 0; i < 3; i++) {
- for (int j = 0; j < 3; j++) {
- mat.setElement(i, j, vec1.v()[i] * vec2.v()[j]);
- }
- }
+
+ private static Matrix v2m(Hep3Vector v) {
+ BasicMatrix mat = new BasicMatrix(3, 1);
+ mat.setElement(0, 0, v.x());
+ mat.setElement(1, 0, v.y());
+ mat.setElement(2, 0, v.z());
return mat;
}
-
- /**
- * Form the outer product matrix vec * vec^T from the vector vec.
- * @param vec vector to be used in creating outer product
- * @return outer product
- */
- private static BasicMatrix Vec2Matrix(Hep3Vector vec) {
- return Vec2Matrix(vec, vec);
- }
-
-
- /**
- * Matrix addition method that is for some reason missing from matrix
- * package. In this method, m2 = m1 + m2.
- * @param m1 first matrix
- * @param m2 second matrix, containing sum when method returns
- */
- private static void MatrixAdd(Matrix m1, MutableMatrix m2) {
- // Add the contents of matrix m1 to matrix m2
- if (m1.getNColumns() != m2.getNColumns()) throw new RuntimeException("HelicalTrackCross matrix dimension mismatch!");
- if (m1.getNRows() != m2.getNRows()) throw new RuntimeException("HelicalTrackCross matrix dimension mismatch");
- for (int i = 0; i < m1.getNRows(); i++) {
- for (int j = 0; j < m1.getNColumns(); j++) {
- m2.setElement(i, j, m1.e(i, j) + m2.e(i, j));
- }
- }
- return;
- }
-
- /**
- * Returns the transpose of the matrix (again, inexplicably not handled by
- * the matrix package for non-square matrices).
- * @param m matrix to be transposed
- * @return transposed matrix
- */
- private static Matrix MatrixTranspose(Matrix m) {
- MutableMatrix mt = new BasicMatrix(m.getNColumns(), m.getNRows());
- for (int i = 0; i < m.getNRows(); i++) {
- for (int j = 0; j < m.getNColumns(); j++) {
- mt.setElement(j, i, m.e(i, j));
- }
- }
- return mt;
- }
-
-}
+}
\ No newline at end of file