1 added + 4 removed + 12 modified, total 17 files
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.2 -r1.3
--- EMap.java 15 Jul 2006 09:12:45 -0000 1.2
+++ EMap.java 7 Aug 2006 21:30:27 -0000 1.3
@@ -1,5 +1,5 @@
/**
- * @version $Id: EMap.java,v 1.2 2006/07/15 09:12:45 jstrube Exp $
+ * @version $Id: EMap.java,v 1.3 2006/08/07 21:30:27 jstrube Exp $
*/
package org.lcsim.contrib.JanStrube.tracking;
@@ -25,10 +25,10 @@
void setPt(double p_t) {
pt = p_t;
}
- double[] getValues() {
+ public double[] getValues() {
return values;
}
- double getPt() {
+ public double getPt() {
return pt;
}
}
\ No newline at end of file
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.4 -r1.5
--- HelixSwimmer.java 18 Jul 2006 03:26:50 -0000 1.4
+++ HelixSwimmer.java 7 Aug 2006 21:30:27 -0000 1.5
@@ -23,7 +23,7 @@
* For more info on swimming see <a href="doc-files/transport.pdf">this paper</a> by Paul Avery.
*
* @author jstrube
- * @version $Id: HelixSwimmer.java,v 1.4 2006/07/18 03:26:50 jstrube Exp $
+ * @version $Id: HelixSwimmer.java,v 1.5 2006/08/07 21:30:27 jstrube Exp $
*/
public class HelixSwimmer {
public class SpatialParameters {
@@ -98,10 +98,39 @@
// System.out.println("Origin is " + origin);
}
- public SpacePoint getPointAtDistance(double alpha) {
+ // FIXME only needed until merger of the two track interfaces is completed
+ /**
+ * @deprecated
+ * @param t
+ */
+ @Deprecated
+ public void setTrack(org.lcsim.event.Track t) {
+ double r = 1/t.getTrackParameter(omega.ordinal());
+ double phi = t.getTrackParameter(phi0.ordinal());
+ double lambda = Math.atan(t.getTrackParameter(tanLambda.ordinal()));
+ double d_0 = t.getTrackParameter(d0.ordinal());
+ double z_0 = t.getTrackParameter(z0.ordinal());
+
+ double[] ref = t.getReferencePoint();
+ double x = ref[0] - d_0 * sin(phi);
+ double y = ref[1] + d_0 * cos(phi);
+ double z = ref[2] + z_0;
+ SpacePoint origin = new CartesianPoint(x, y, z);
+ _trajectory = new Helix(origin, r, phi, lambda);
+
+ double pt = r*t.getCharge()*field;
+ spatialParms.px = pt*cos(phi);
+ spatialParms.py = pt*sin(phi);
+ spatialParms.pz = pt*t.getTrackParameter(tanLambda.ordinal());
+ spatialParms.charge = t.getCharge();
+// System.out.println(spatialParms);
+// System.out.println("Origin is " + origin);
+ }
+
+ public SpacePoint getPointAtLength(double alpha) {
if (_trajectory == null)
throw new RuntimeException("Trajectory not set");
- return _trajectory.getPointAtDistance(alpha);
+ return _trajectory.getPointAtLength(alpha);
}
public double getDistanceToRadius(double r) {
@@ -135,15 +164,21 @@
return _trajectory.getTrackLengthToPoint(point);
}
+ // FIXME legacy implementation
+ @Deprecated
+ public double getTrackLengthToPoint(double[] x) {
+ return _trajectory.getTrackLengthToPoint(new CartesianPoint(x));
+ }
+
/**
* Returns the momentum on a point on the track at a distance from the origin
*
* @param alpha The 3D distance from the origin along the track
* @return The components of the momentum in a SpacePoint
*/
- public SpacePoint getMomentumAtDistance(double alpha) {
+ public SpacePoint getMomentumAtLength(double alpha) {
// the trajectory can only return the unit direction of the momentum
- SpacePoint unitDirection = _trajectory.getMomentumAtLength(alpha);
+ SpacePoint unitDirection = _trajectory.getUnitTangentAtLength(alpha);
return VectorArithmetic.multiply(unitDirection, sqrt(spatialParms.px*spatialParms.px + spatialParms.py*spatialParms.py));
}
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.1 -r1.2
--- Line.java 17 Jul 2006 10:12:39 -0000 1.1
+++ Line.java 7 Aug 2006 21:30:27 -0000 1.2
@@ -11,7 +11,7 @@
* A straight line
*
* @author tonyj
- * @version $Id: Line.java,v 1.1 2006/07/17 10:12:39 jstrube Exp $
+ * @version $Id: Line.java,v 1.2 2006/08/07 21:30:27 jstrube Exp $
*/
public class Line implements Trajectory {
/**
@@ -325,7 +325,7 @@
/**
* Gets a point after traveling distance alpha from the trajectory's origin point along the trajectory
*/
- public SpacePoint getPointAtDistance(double alpha) {
+ public SpacePoint getPointAtLength(double alpha) {
double z = origin.z() + alpha * sinLambda;
double dr = alpha * cosLambda;
double x = origin.x() + dr * cosPhi;
@@ -350,7 +350,7 @@
* @param alpha is ignored
* @return The unit direction of the line, since there is now geometric way to obtain the momentum along a straight line.
*/
- public SpacePoint getMomentumAtLength(double alpha) {
+ public SpacePoint getUnitTangentAtLength(double alpha) {
SpacePoint lineDirection = new CartesianPoint(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda);
return lineDirection;
}
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.4 -r1.5
--- NewTrack.java 3 Aug 2006 10:35:50 -0000 1.4
+++ NewTrack.java 7 Aug 2006 21:30:27 -0000 1.5
@@ -1,12 +1,27 @@
/**
- * @version $Id: NewTrack.java,v 1.4 2006/08/03 10:35:50 jstrube Exp $
+ * @version $Id: NewTrack.java,v 1.5 2006/08/07 21:30:27 jstrube Exp $
*/
package org.lcsim.contrib.JanStrube.tracking;
+import org.lcsim.spacegeom.CartesianPoint;
import org.lcsim.spacegeom.SpacePoint;
import Jama.Matrix;
+import static org.lcsim.contrib.JanStrube.tracking.Track.ParameterName.d0;
+import static org.lcsim.contrib.JanStrube.tracking.Track.ParameterName.phi0;
+import static org.lcsim.contrib.JanStrube.tracking.Track.ParameterName.omega;
+import static org.lcsim.contrib.JanStrube.tracking.Track.ParameterName.z0;
+import static org.lcsim.contrib.JanStrube.tracking.Track.ParameterName.tanLambda;
+
+import static org.lcsim.contrib.JanStrube.tracking.Constants.fieldConversion;
+
+import static java.lang.Math.sin;
+import static java.lang.Math.cos;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.atan2;
+
+
/**
* @author jstrube
* The class to store the measurement information of the track of a charged particle in a magnetic field.
@@ -32,21 +47,57 @@
}
public static SpacePoint Parameters2Position(EMap parameters) {
- return Parameters2Position(parameters, new SpacePoint());
+ final SpacePoint zero = new SpacePoint();
+ return Parameters2Position(parameters, zero);
}
-
+
public static SpacePoint Parameters2Momentum(EMap parameters) {
- return Parameters2Momentum(parameters, new SpacePoint());
- }
-
- public static SpacePoint Parameters2Momentum(EMap parameters, SpacePoint refPoint) {
- return new SpacePoint();
+ double pt = parameters.pt;
+
+ return new CartesianPoint(
+ pt * sin(parameters.get(phi0))
+ , pt * cos(parameters.get(phi0))
+ , pt * cos(parameters.get(tanLambda))
+ );
}
public static SpacePoint Parameters2Position(EMap parameters, SpacePoint refPoint) {
- return new SpacePoint();
+ double d_0 = parameters.get(d0);
+ double z_0 = parameters.get(z0);
+ double phi_0 = parameters.get(phi0);
+
+ double x = refPoint.x() + d_0*sin(phi_0);
+ double y = refPoint.y() - d_0*sin(phi_0);
+ double z = refPoint.z() + z_0;
+
+ return new CartesianPoint(x, y, z);
+ }
+
+ public static EMap SpaceMomentum2Parameters(SpacePoint pos, SpacePoint mom, SpacePoint ref, int charge, double Bz) {
+ EMap result = new EMap();
+ double aqBz = charge*Bz*fieldConversion;
+ double x = pos.x() - ref.x();
+ double y = pos.y() - ref.y();
+ double z = pos.z() - ref.z();
+ double px0 = mom.x() - aqBz*y;
+ double py0 = mom.y() + aqBz*x;
+ double pt0 = sqrt(px0*px0 + py0*py0);
+
+ double phi_0 = atan2(py0, px0);
+ double phi = mom.phi();
+ double pt = mom.rxy();
+ double l = (phi_0 - phi)*pt/aqBz;
+
+ result.set(d0, (pt-pt0)/aqBz);
+ result.set(phi0, phi_0);
+ result.set(omega, aqBz/pt);
+ result.set(z0, z - l*mom.z()/pt);
+ result.set(tanLambda, mom.z()/pt);
+ result.setPt(pt);
+ return result;
}
+
/* (non-Javadoc)
* @see org.lcsim.contrib.JanStrube.tracking.Track#getParameter(org.lcsim.contrib.JanStrube.tracking.NewTrack.ParameterName)
*/
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.6 -r1.7
--- NewFastMCTrackFactory.java 19 Jul 2006 23:14:01 -0000 1.6
+++ NewFastMCTrackFactory.java 7 Aug 2006 21:30:27 -0000 1.7
@@ -1,5 +1,5 @@
/**
- * @version $Id: NewFastMCTrackFactory.java,v 1.6 2006/07/19 23:14:01 jstrube Exp $
+ * @version $Id: NewFastMCTrackFactory.java,v 1.7 2006/08/07 21:30:27 jstrube Exp $
*/
package org.lcsim.contrib.JanStrube.tracking;
@@ -130,8 +130,8 @@
Track getTrack(SpacePoint momentum, SpacePoint location, SpacePoint referencePoint, int charge, Random random, boolean shouldISmear) {
_swimmer.setTrack(momentum, location, charge);
double alpha = _swimmer.getTrackLengthToPoint(referencePoint);
- SpacePoint poca = _swimmer.getPointAtDistance(alpha);
- SpacePoint momentumAtPoca = _swimmer.getMomentumAtDistance(alpha);
+ SpacePoint poca = _swimmer.getPointAtLength(alpha);
+ SpacePoint momentumAtPoca = _swimmer.getMomentumAtLength(alpha);
EMap parameters = new EMap();
parameters.pt = momentumAtPoca.rxy();
parameters.set(phi0, atan2(momentumAtPoca.y(), momentumAtPoca.x()));
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.1 -r1.2
--- Trajectory.java 17 Jul 2006 10:12:39 -0000 1.1
+++ Trajectory.java 7 Aug 2006 21:30:27 -0000 1.2
@@ -7,14 +7,14 @@
/**
* A particle trajectory (either a Helix or a Line)
* @author tonyj
- * @version $Id: Trajectory.java,v 1.1 2006/07/17 10:12:39 jstrube Exp $
+ * @version $Id: Trajectory.java,v 1.2 2006/08/07 21:30:27 jstrube Exp $
*/
public interface Trajectory
{
/**
* Gets a point after traveling distance alpha from the origin along the trajectory
*/
- SpacePoint getPointAtDistance(double alpha);
+ SpacePoint getPointAtLength(double alpha);
/**
* Calculates the distance at which the trajectory first reaches radius R.
* Returns Double.NaN if the trajectory does not intercept the cylinder.
@@ -43,5 +43,5 @@
* @param alpha The length along the trajectory from the origin
* @return The momentum at the given distance
*/
- public SpacePoint getMomentumAtLength(double alpha);
+ public SpacePoint getUnitTangentAtLength(double alpha);
}
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.3 -r1.4
--- Helix.java 19 Jul 2006 23:14:01 -0000 1.3
+++ Helix.java 7 Aug 2006 21:30:27 -0000 1.4
@@ -22,7 +22,7 @@
* For more info on swimming see <a href="doc-files/transport.pdf">this paper</a> by Paul Avery.
*
* @author tonyj
- * @version $Id: Helix.java,v 1.3 2006/07/19 23:14:01 jstrube Exp $
+ * @version $Id: Helix.java,v 1.4 2006/08/07 21:30:27 jstrube Exp $
*/
public class Helix implements Trajectory {
private Hep3Vector origin;
@@ -107,7 +107,7 @@
* @param alpha The distance along the trajectory
* @return The unit vector of the momentum
*/
- public SpacePoint getMomentumAtLength(double alpha) {
+ public SpacePoint getUnitTangentAtLength(double alpha) {
double angle = phi + alpha * rho;
return new CartesianPoint(cos(angle), sin(angle), sinLambda / cosLambda);
}
@@ -116,7 +116,7 @@
* returns a point on the Helix at a distance alpha from the origin along the trajectory. alpha ==
* 2*PI*radius/cos(lambda) is one rotation in the x-y plane
*/
- public SpacePoint getPointAtDistance(double alpha) {
+ public SpacePoint getPointAtLength(double alpha) {
double angle = phi + alpha * rho;
double x = xCenter - radius * sin(angle);
double y = yCenter + radius * cos(angle);
@@ -163,14 +163,15 @@
return numerator / denominator;
}
- // Returns the "momentum" at the length s from the starting point.
- // This uses the definition in http://www.phys.ufl.edu/~avery/fitting/transport.pdf
- public SpacePoint getTangentAtDistance(double alpha) {
- double p0x = px * cos(rho * alpha) - py * sin(rho * alpha);
- double p0y = py * cos(rho * alpha) + px * sin(rho * alpha);
- double p0z = pz * cos(rho * alpha) + pz * (1 - cos(rho * alpha));
- return new CartesianPoint(p0x, p0y, p0z);
- }
+ // FIXME this method seems obsolete -- remove it.
+// // Returns the "momentum" at the length s from the starting point.
+// // This uses the definition in http://www.phys.ufl.edu/~avery/fitting/transport.pdf
+// public SpacePoint getTangentAtLength(double alpha) {
+// double p0x = px * cos(rho * alpha) - py * sin(rho * alpha);
+// double p0y = py * cos(rho * alpha) + px * sin(rho * alpha);
+// double p0z = pz * cos(rho * alpha) + pz * (1 - cos(rho * alpha));
+// return new CartesianPoint(p0x, p0y, p0z);
+// }
/**
* Swims the helix along its trajectory to the point of closest approach to the given SpacePoint. For more info on
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -u -r1.5 -r1.6
--- Fitter.java 4 Jul 2006 01:41:48 -0000 1.5
+++ Fitter.java 7 Aug 2006 21:30:30 -0000 1.6
@@ -1,47 +1,42 @@
package org.lcsim.contrib.JanStrube.vtxFitter;
/**
- * @version $Id: Fitter.java,v 1.5 2006/07/04 01:41:48 jstrube Exp $
+ * @version $Id: Fitter.java,v 1.6 2006/08/07 21:30:30 jstrube Exp $
*/
-import org.lcsim.event.ReconstructedParticle;
-import org.lcsim.event.Track;
+import static java.lang.Math.atan2;
+import static java.lang.Math.sqrt;
+import static org.lcsim.contrib.JanStrube.tracking.Constants.fieldConversion;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.subtract;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.add;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.dot;
+import static org.lcsim.contrib.JanStrube.tracking.NewTrack.SpaceMomentum2Parameters;
+
+import org.lcsim.contrib.JanStrube.tracking.HelixSwimmer;
+import org.lcsim.contrib.JanStrube.tracking.Track;
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.SpacePoint;
import Jama.Matrix;
public class Fitter {
- ReconstructedParticle particle;
+ // The Vertex object
+ Vertex particle;
+
+ Vertex unfittedParticle;
// convergence criterion
double deltaChi2 = 1e-4;
+ double B_z;
// position
- double[] x_prev;
+ Matrix x_prev;
// covariance for x
- Matrix C_prev;
- // momentum
- double[] q_prev;
- // covariance for momentum
- Matrix D_prev;
- //covariance for x and q
- Matrix E_prev;
- // covariance for track measurement
- Matrix V_prev;
- // inverse of covariance matrix
- Matrix G_prev;
-
- // linearity factors of the measeurement eq.
- // in x
- Matrix A_k, A_prev, A;
- // and q
- Matrix B_k, B_prev, B;
-
- // residual
- double[] c_k0;
+ Matrix C_prev;
+
+ HelixSwimmer _swimmer;
// chi2
double chi2_prev;
-
- double[] r_k;
-
+
/**
* Definitions:
* x_k = estimate of the vertex position after using the information of k tracks
@@ -57,56 +52,116 @@
* Gk = Vk^-1 weightmatrix of particle k
* @param p
*/
- public Fitter(ReconstructedParticle p) {
+ public Fitter(Vertex p, double bz) {
this(p
- , new double[3]
- , new double[][] {{10000, 0, 0}, {0, 10000, 0}, {0, 0, 10000}});
+ , bz
+ , new Matrix(new double[][] {{10000, 0, 0}, {0, 10000, 0}, {0, 0, 10000}})
+ );
}
- public Fitter(ReconstructedParticle p, double[] x0, double[][] C0) {
- particle = p;
- x_prev = x0;
- C_prev = new Matrix(C0);
+ public Fitter(Vertex p, double bz, Matrix c0) {
+ unfittedParticle = p;
+ B_z = bz;
+ C_prev = c0;
+ x_prev = new Matrix(p.origin().getCartesianArray(), 3);
+ _swimmer = new HelixSwimmer(bz);
}
- public void fit() {
- for (Track t : particle.getTracks()) {
+ public Vertex fit() {
+ particle = new Vertex();
+ particle.setOrigin(unfittedParticle.origin());
+
+ for (Track t : unfittedParticle.getTracks()) {
filter(t);
}
for (Track t : particle.getTracks()) {
smoothe(t);
}
+ return particle;
}
+
+ // The problem is that double[] is a column vector...
+ // Turn everything into a Matrix
private double filter(Track t) {
- Matrix G_k = new Matrix(t.getErrorMatrix()).inverse();
- double[] m_k = t.getTrackParameters();
-
- Matrix W_k = (B_k.transpose().times(G_k.times(B_k))).uminus();
- Matrix G_B = G_k.minus((G_k.times(B_k.times(W_k.times(B_k.transpose().times(G_k))))));
- Matrix C_k = C_prev.inverse().plus(A_k.transpose().times(G_B.times(A_k))).inverse();
-
- double[] x_k = C_k.times(add(C_prev.inverse().times(x_prev)
- , A_k.transpose().times(G_B).times(subtract(m_k, c_k0))));
+ // First, get the Matrices at the point of linearization x0, p0
+ // a good choice for x0 is the previous best estimate and
+ // p0 is the momentum of the track at the POCA to x0
+ _swimmer.setTrack(t);
+ // at the moment, particle.origin is the best estimate so far
+ double alpha = _swimmer.getTrackLengthToPoint(particle.origin());
+ SpacePoint p = _swimmer.getMomentumAtLength(alpha);
+ // x is the POCA of the current track to the best estimate
+
+ SpacePoint x = _swimmer.getPointAtLength(alpha);
+
+ Matrix A_k = getSpatialDerivativeMatrix(particle.origin(), p, t.getCharge(), B_z);
+ Matrix B_k = getMomentumDerivativeMatrix(particle.origin(), p, t.getCharge(), B_z);
+ Matrix q0 = new Matrix(p.getCartesianArray(), 3);
+ // turn the measurement into a Matrix object
+ Matrix h = new Matrix(SpaceMomentum2Parameters(
+ particle.origin(), p, particle.origin(), t.getCharge(), B_z).getValues()
+ , 5);
+ // the measurement of the track is taken at the POCA to the best estimate
+ Matrix m_k = new Matrix(SpaceMomentum2Parameters(
+ x, p, particle.origin(), t.getCharge(), B_z).getValues()
+ , 5);
+ Matrix c_k0 = h.minus(A_k.times(x_prev).plus(B_k.times(q0)));
+
+ // Now for the real filtering equations
+ Matrix G_k = t.getErrorMatrix().inverse();
+ Matrix W_k = (B_k.transpose().times(G_k.times(B_k))).inverse();
+ Matrix G_kB = G_k.minus((G_k.times(B_k.times(W_k.times(B_k.transpose().times(G_k))))));
+ // cov(x)
+ Matrix C_k = C_prev.inverse().plus(A_k.transpose().times(G_kB.times(A_k))).inverse();
+ Matrix x_k = C_k.times(C_prev.inverse().times(x_prev).plus(
+ A_k.transpose().times(G_kB).times(m_k.minus(c_k0))));
+ Matrix q_k = W_k.times(B_k.transpose().times(G_k.times(m_k.minus(c_k0.plus(A_k.times(x_k))))));
+
+ // cov(q)
+ Matrix D_k = W_k.plus(W_k.times(B_k.transpose().times(G_k.times(A_k.times(C_k.times(A_k.transpose().times(G_k.times(B_k.times(W_k)))))))));
+ // cov(x, q)
+ Matrix E_k = C_k.times(A_k.transpose().times(G_k.times(B_k.times(W_k)))).uminus();
+
+ // now the chi2
+ Matrix r_k = h.minus(c_k0.plus(A_k.times(x_k).plus(B_k.times(q_k))));
+ double chi2_kf = r_k.transpose().times(G_k.times(r_k)).get(0, 0)
+ + x_k.minus(x_prev).transpose().times(C_prev.inverse().times(x_k.minus(x_prev))).get(0, 0);
- double[] q_k = W_k.times(B_k.transpose().times(G_k)).times(subtract(m_k, add(c_k0, A_k.times(x_k))));
- Matrix D_k = W_k.plus(W_k.times(B_k.transpose().plus(G_k.times(A_k.times(C_k.times(A_k.transpose().times(G_k.times(B_k.times(W_k)))))))));
- Matrix E_k = C_k.times(A_k.transpose().times(G_k.times(B_k.times(W_k)))).uminus();
-
-
- double chi2 = chi2_prev + dot(r_k, G_k.times(r_k)) + dot(subtract(x_k, x_prev), C_prev.inverse().times(subtract(x_k, x_prev)));
+
+ double chi2 = chi2_prev + chi2_kf;
- r_k = subtract(m_k, add(c_k0, add(A_k.times(x_k), B_k.times(q_k))));
+ x_prev = x_k;
C_prev = C_k;
+ particle.addTrack(t, chi2_kf);
+ particle._origin = new CartesianPoint(x_k.getRowPackedCopy());
+ chi2_prev = chi2;
return chi2;
}
private void smoothe(Track t) {
- Matrix G_k = new Matrix(t.getErrorMatrix()).inverse();
- Matrix W_k = (B_k.transpose().times(G_k.times(B_k))).uminus();
- double[] m_k = t.getTrackParameters();
- double[] q_kN = W_k.times(B_k.transpose().times(G_k)).times(subtract(m_k, add(c_k0, A_k.times(x_prev))));
+ _swimmer.setTrack(t);
+ // at the moment, particle.origin is the best estimate so far
+ double alpha = _swimmer.getTrackLengthToPoint(particle.origin());
+ SpacePoint p = _swimmer.getMomentumAtLength(alpha);
+ // x is the POCA of the current track to the best estimate
+ Matrix q0 = new Matrix(p.getCartesianArray(), 3);
+
+ SpacePoint x = _swimmer.getPointAtLength(alpha);
+
+ Matrix A_k = getSpatialDerivativeMatrix(particle.origin(), p, t.getCharge(), B_z);
+ Matrix B_k = getMomentumDerivativeMatrix(particle.origin(), p, t.getCharge(), B_z);
+
+ Matrix G_k = t.getErrorMatrix().inverse();
+ Matrix W_k = B_k.transpose().times(G_k.times(B_k)).uminus();
+ // the measurement of the track is taken at the POCA to the best estimate
+ Matrix m_k = new Matrix(SpaceMomentum2Parameters(
+ x, p, particle.origin(), t.getCharge(), B_z).getValues()
+ , 5);
+ Matrix c_k0 = m_k.minus(A_k.times(x_prev).plus(B_k.times(q0)));
+
+ Matrix q_kN = W_k.times(B_k.transpose().times(G_k)).times(m_k.minus(c_k0.plus(A_k.times(x_prev))));
Matrix D_kN = W_k.plus(W_k.times(B_k.transpose().times(G_k.times(A_k.times(C_prev.times(A_k.transpose().times(G_k.times(B_k.times(W_k)))))))));
Matrix E_kN = C_prev.times(A_k.transpose().times(G_k.times(B_k.times(W_k)))).uminus();
}
@@ -130,26 +185,93 @@
// return a.times(b);
// }
- static double[] subtract(double[] a, double[] b) {
- if (a.length != b.length)
- throw new IllegalArgumentException("dimensions do not match");
- double[] result = new double[a.length];
- for (int i=0; i<a.length; i++) {
- result[i] = a[i] - b[i];
- }
- return result;
- }
-
- static double[] add(double[] a, double[] b) {
- if (a.length != b.length)
- throw new IllegalArgumentException("dimensions do not match");
- double[] result = new double[a.length];
- for (int i=0; i<a.length; i++) {
- result[i] = a[i] + b[i];
- }
- return result;
- }
+ public static Matrix getMomentumDerivativeMatrix(SpacePoint x, SpacePoint p, int charge, double Bz) {
+ double field = Bz*fieldConversion;
+ double pt = p.rxy();
+ double px0 = p.x() - charge*field;
+ double py0 = p.y() + charge*field;
+ double pt0 = sqrt(px0*px0 + py0*py0);
+ double phi = atan2(p.y(), p.x());
+ double phi0 = atan2(py0, px0);
+ double l = (phi-phi0)*pt/(charge*field);
+
+ double[] dTanLambda_dp = new double[] {
+ -p.z()*p.x() / (pt*pt*pt)
+ , -p.z()*p.y() / (pt*pt*pt)
+ , 1/pt
+ };
+ double[] dPhi_dp = new double[] {
+ -py0 / (pt0*pt0)
+ , px0 / (pt0*pt0)
+ , 0
+ };
+ double[] dOmega_dp = new double[] {
+ -charge*field*p.x() / (pt*pt*pt)
+ , -charge*field*p.y() / (pt*pt*pt)
+ , 0
+ };
+ double[] dD0_dp = new double[] {
+ -py0/pt0
+ , px0/pt0
+ , 0
+ };
+ double[] dZ_dp = new double[] {
+ -p.z() / (field*charge) * (p.y()/(pt*pt) - py0/(pt0*pt0))
+ , -p.z() / (field*charge) * (p.x()/(pt*pt) - px0/(pt0*pt0))
+ , -l/pt
+ };
+ return new Matrix(new double[][] {
+ dD0_dp, dPhi_dp, dOmega_dp, dZ_dp, dTanLambda_dp
+ });
+ }
+ /**
+ * Calculates the derivatives of the measurement vector with respect to cartesian spatial coordinates.
+ * @param x
+ * @param p
+ * @param charge
+ * @return the derivatives Matrix. Assumes an ordering of ??? in the measurement vector
+ */
+ public static Matrix getSpatialDerivativeMatrix(SpacePoint x, SpacePoint p, int charge, double Bz) {
+ double field = Bz*fieldConversion;
+
+ double px0 = p.x() - charge*field;
+ double py0 = p.y() + charge*field;
+ double pt0 = sqrt(px0*px0 + py0*py0);
+
+ // be very explicit about the derivatives
+ // dr = (dx, dy, dz)^T
+ double[] dTanLambda_dr = new double[] {
+ 0, 0, 0
+ };
+
+ double[] dPhi_dr = new double[] {
+ px0*field*charge / (pt0*pt0)
+ , py0*field*charge / (pt0*pt0)
+ , 0
+ };
+
+ double[] dOmega_dr = new double[] {
+ 0, 0, 0
+ };
+
+ double[] dD0_dr = new double[] {
+ -py0 / pt0
+ , px0 / pt0
+ , 0
+ };
+
+ double[] dZ_dr = new double[] {
+ -p.z() * px0 / (pt0*pt0)
+ , -p.z() * py0 / (pt0*pt0)
+ , 1
+ };
+
+ return new Matrix(new double[][] {
+ dD0_dr, dPhi_dr, dOmega_dr, dZ_dr, dTanLambda_dr
+ });
+ }
+
// for notational convenience
static Matrix add(Matrix a, Matrix b) {
return a.plus(b);
@@ -164,15 +286,6 @@
static Matrix chargedMomentumLinearity(Track t) {
return new Matrix(0, 0);
}
-
- static double dot(double[] a, double[] b) {
- if (a.length != b.length)
- throw new IllegalArgumentException("dimensions don't match");
- double result = 0;
- for (int i=0; i<a.length; i++)
- result += a[i]*b[i];
- return result;
- }
static Matrix calcSpatialLinearityMatrix(Track t) {
return new Matrix(5, 5);
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -u -r1.1 -r1.2
--- Vertex.java 28 Jun 2006 01:54:49 -0000 1.1
+++ Vertex.java 7 Aug 2006 21:30:30 -0000 1.2
@@ -1,32 +1,73 @@
package org.lcsim.contrib.JanStrube.vtxFitter;
-import org.lcsim.event.base.BaseReconstructedParticle;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import org.lcsim.contrib.JanStrube.tracking.Track;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.ReconstructedParticle;
import org.lcsim.spacegeom.CartesianPoint;
import org.lcsim.spacegeom.SpacePoint;
-public class Vertex extends BaseReconstructedParticle {
- private SpacePoint origin;
- private SpacePoint endpoint;
-
- public Vertex(double mass) {
- super(mass);
- origin = new CartesianPoint(0, 0, 0);
- endpoint = new CartesianPoint(0, 0, 0);
+import Jama.Matrix;
+
+public class Vertex {
+ SpacePoint _origin;
+ private SpacePoint _originError;
+ private Map<Track, Double> _tracks;
+
+ double _chi2;
+ private double _ndf;
+ private Matrix _covMatrix;
+
+ public Vertex() {
+ _origin = new SpacePoint();
+ _originError = new SpacePoint();
+ _tracks = new HashMap<Track, Double>();
}
public void setOrigin(SpacePoint x) {
- origin = x;
- }
-
- public void setEndPoint(SpacePoint x) {
- endpoint = x;
- }
-
- public SpacePoint endpoint() {
- return endpoint;
+ _origin = x;
}
public SpacePoint origin() {
- return origin;
+ return _origin;
}
+
+ public void addTrack(Track t, double chi2) {
+ _tracks.put(t, chi2);
+ }
+
+ public void setTracks(Map<Track, Double> tracks) {
+ _tracks = tracks;
+ }
+
+ public void setTracks(List<Track> tracks) {
+ for (Track t : tracks) {
+ _tracks.put(t, 0.);
+ }
+ }
+
+ public double getChi2() {
+ return _chi2;
+ }
+
+ public void setChi2(double chi2) {
+ _chi2 = chi2;
+ }
+
+ public double getNdf() {
+ return _ndf;
+ }
+
+ public void setNdf(double ndf) {
+ _ndf = ndf;
+ }
+
+ public Set<Track> getTracks() {
+ return _tracks.keySet();
+ }
}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -u -r1.5 -r1.6
--- vertexing.lyx 3 Aug 2006 10:35:53 -0000 1.5
+++ vertexing.lyx 7 Aug 2006 21:30:30 -0000 1.6
@@ -261,7 +261,7 @@
\begin_layout Standard
\begin_inset Formula \begin{align*}
-pt & =\sqrt{p_{x}^{2}+p_{y}^{2}}\\
+p_{t} & =\sqrt{p_{x}^{2}+p_{y}^{2}}\\
\omega & =aq\frac{B_{z}}{p_{t}}\\
\phi_{0} & =\text{atan2}(py,px)\\
\tan(\lambda) & =\frac{p_{z}}{p_{t}}\\
@@ -495,7 +495,15 @@
\begin_layout Standard
With these conventions established, the derivatives wrt.
- the components of space and momentum are:
+ the components of space and momentum are (substituting
+\begin_inset Formula $\omega=aqB_{z}/p_{t}$
+\end_inset
+
+ to avoid an explicit dependence on
+\begin_inset Formula $B_{z}$
+\end_inset
+
+):
\end_layout
\begin_layout Standard
@@ -504,24 +512,24 @@
\frac{\dd\tan(\lambda)}{\dd p_{x}} & =-\frac{p_{z}p_{x}}{p_{t}^{3}}\\
\frac{\dd\tan(\lambda)}{\dd p_{y}} & =-\frac{p_{z}p_{y}}{p_{t}^{3}}\\
\frac{\dd\tan(\lambda)}{\dd p_{z}} & =\frac{1}{p_{t}}\\
-\frac{\dd\phi_{0}}{\dd x} & =\frac{\tilde{p}_{x}aqB_{z}}{\tilde{p}_{t}^{2}}\\
-\frac{\dd\phi_{0}}{\dd y} & =\frac{\tilde{p}_{y}aqB_{z}}{\tilde{p}_{t}^{2}}\\
+\frac{\dd\phi_{0}}{\dd x} & =\frac{\tilde{p}_{x}aqB_{z}}{\tilde{p}_{t}^{2}}=\frac{\tilde{p}_{x}\omega p_{t}}{\tilde{p}_{t}^{2}}\\
+\frac{\dd\phi_{0}}{\dd y} & =\frac{\tilde{p}_{y}aqB_{z}}{\tilde{p}_{t}^{2}}=\frac{\tilde{p}_{y}\omega p_{t}}{\tilde{p}_{t}^{2}}\\
\frac{\dd\phi_{0}}{\dd z} & =\frac{\dd\phi_{0}}{\dd p_{z}}=0\\
\frac{\dd\phi_{0}}{\dd p_{x}} & =-\frac{\tilde{p}_{y}}{\tilde{p}_{t}^{2}}\\
\frac{\dd\phi_{0}}{\dd p_{y}} & =\frac{\tilde{p}_{x}}{\tilde{p}_{t}^{2}}\\
\frac{\dd\omega}{\dd x} & =\frac{\dd\omega}{\dd y}=\frac{\dd\omega}{\dd z}=\frac{\dd\omega}{\dd p_{z}}=0\\
-\frac{\dd\omega}{\dd p_{x}} & =-\frac{aqB_{z}p_{x}}{p_{t}^{3}}\\
-\frac{\dd\omega}{\dd p_{y}} & =-\frac{aqB_{z}p_{y}}{p_{t}^{3}}\\
+\frac{\dd\omega}{\dd p_{x}} & =-\frac{aqB_{z}p_{x}}{p_{t}^{3}}=-\frac{\omega p_{x}}{p_{t}^{2}}\\
+\frac{\dd\omega}{\dd p_{y}} & =-\frac{aqB_{z}p_{y}}{p_{t}^{3}}=-\frac{\omega p_{y}}{p_{t}^{2}}\\
\frac{\dd d_{0}}{\dd x} & =-\frac{\tilde{p}_{y}}{\tilde{p}_{t}}\\
\frac{\dd d_{0}}{\dd y} & =\frac{\tilde{p}_{x}}{\tilde{p}_{t}}\\
\frac{\dd d_{0}}{\dd z} & =\frac{\dd d_{0}}{\dd p_{z}}=0\\
-\frac{\dd d_{0}}{\dd p_{x}} & =\frac{1}{aqB_{z}}(\frac{p_{x}}{p_{t}}-\frac{\tilde{p}_{x}}{\tilde{p}_{t}})\\
-\frac{\dd d_{0}}{\dd p_{y}} & =\frac{1}{aqB_{z}}(\frac{p_{y}}{p_{t}}-\frac{\tilde{p}_{y}}{\tilde{p}_{t}})\\
+\frac{\dd d_{0}}{\dd p_{x}} & =\frac{1}{aqB_{z}}(\frac{p_{x}}{p_{t}}-\frac{\tilde{p}_{x}}{\tilde{p}_{t}})=\frac{1}{\omega p_{t}}(\frac{p_{x}}{p_{t}}-\frac{\tilde{p}_{x}}{\tilde{p}_{t}})\\
+\frac{\dd d_{0}}{\dd p_{y}} & =\frac{1}{aqB_{z}}(\frac{p_{y}}{p_{t}}-\frac{\tilde{p}_{y}}{\tilde{p}_{t}})=\frac{1}{\omega p_{t}}(\frac{p_{y}}{p_{t}}-\frac{\tilde{p}_{y}}{\tilde{p}_{t}})\\
\frac{\dd z_{0}}{\dd x} & =-\frac{p_{z}\tilde{p}_{x}}{\tilde{p}_{t}^{2}}\\
\frac{\dd z_{0}}{\dd y} & =-\frac{p_{z}\tilde{p}_{y}}{\tilde{p}_{t}^{2}}\\
\frac{\dd z_{0}}{\dd z} & =1\\
-\frac{\dd z_{0}}{\dd p_{x}} & =-\frac{p_{z}}{aqB_{z}}(\frac{p_{y}}{p_{t}^{2}}-\frac{\tilde{p}_{y}}{\tilde{p}_{t}^{2}})\\
-\frac{\dd z_{0}}{\dd p_{y}} & =\frac{p_{z}}{aqB_{z}}(\frac{p_{x}}{p_{t}^{2}}-\frac{\tilde{p}_{x}}{\tilde{p}_{t}^{2}})\\
+\frac{\dd z_{0}}{\dd p_{x}} & =-\frac{p_{z}}{aqB_{z}}(\frac{p_{y}}{p_{t}^{2}}-\frac{\tilde{p}_{y}}{\tilde{p}_{t}^{2}})=-\frac{p_{z}}{\omega p_{t}}(\frac{p_{y}}{p_{t}^{2}}-\frac{\tilde{p}_{y}}{\tilde{p}_{t}^{2}})\\
+\frac{\dd z_{0}}{\dd p_{y}} & =\frac{p_{z}}{aqB_{z}}(\frac{p_{x}}{p_{t}^{2}}-\frac{\tilde{p}_{x}}{\tilde{p}_{t}^{2}})=\frac{p_{z}}{\omega p_{t}}(\frac{p_{x}}{p_{t}^{2}}-\frac{\tilde{p}_{x}}{\tilde{p}_{t}^{2}})\\
\frac{\dd z_{0}}{\dd p_{z}} & =-\frac{l}{p_{t}}\end{align*}
\end_inset
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -N VertexFactory.java
--- VertexFactory.java 18 Jul 2006 03:26:52 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,58 +0,0 @@
-package org.lcsim.contrib.JanStrube.vtxFitter;
-
-import hep.physics.vec.Hep3Vector;
-
-import org.lcsim.contrib.JanStrube.tracking.NewFastMCTrackFactory;
-import org.lcsim.event.ReconstructedParticle;
-import org.lcsim.event.Track;
-import org.lcsim.event.base.BaseReconstructedParticle;
-import org.lcsim.recon.tracking.trfutil.RandomGenerator;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.subtract;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.SpacePoint;
-
-/**
- * Class to create ReconstructedParticle objects consisting
- * of tracks with random (within limits) properties.
- * This isn't very useful at the moment, as the tracks in the vertex are truly random,
- * and not jetlike.
- * @author jstrube
- *
- */
-// TODO This might be interesting to interface with HepRep
-public class VertexFactory extends RandomGenerator {
- SpacePoint momentum;
- int charge;
- /**
- * Creates a ReconstructedParticle at the given position and with the given momentum.
- * The momentum of the daughter tracks is random, but adds up to the given momentum.
- * The charge of the vertex is related to the number of tracks as follows:
- * charge = (-1)^(round(nTracks/2))
- * @param pos the position of the vertex in space
- * @param mom the sum of the momenta of the tracks
- * @param charge the sum of the charges of the tracks
- * @param nTracks the number of tracks
- * @return a new Vertex with the desired properties
- */
- public Vertex makeVertex(SpacePoint pos, SpacePoint mom, int nTracks) {
- if (nTracks < 1)
- throw new IllegalArgumentException("nTracks should be > 1");
- momentum = (SpacePoint) mom.clone();
- Vertex particle = new Vertex(0.);
- particle.setOrigin(pos);
-// NewFastMCTrackFactory fac = new NewFastMCTrackFactory("sidaug05", 5., true);
-// for (int i=0; i<nTracks-1; i++) {
-// int c = (int) Math.pow(-1, Math.round((i+1)/2));
-// double px = flat(-momentum.x(), momentum.x());
-// double py = flat(-momentum.y(), momentum.y());
-// double pz = flat(-momentum.z(), momentum.z());
-// SpacePoint p = new CartesianPoint(px, py, pz);
-// momentum = subtract(momentum, p);
-// Track t = fac.getMCTrack(p, pos, c);
-// particle.addTrack(t);
-// }
-// int c = (int)Math.pow(-1, Math.round((nTracks+1)/2));
-// particle.addTrack(fac.getMCTrack(momentum, pos, c));
- return particle;
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -N VertexTrack.java
--- VertexTrack.java 3 Aug 2006 10:35:53 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,121 +0,0 @@
-/**
- * @version $Id: VertexTrack.java,v 1.5 2006/08/03 10:35:53 jstrube Exp $
- */
-package org.lcsim.contrib.JanStrube.vtxFitter;
-
-import static java.lang.Math.atan2;
-import static java.lang.Math.cos;
-import static java.lang.Math.sin;
-import static java.lang.Math.sqrt;
-import static java.lang.Math.atan;
-import static org.lcsim.contrib.JanStrube.tracking.Constants.fieldConversion;
-
-import org.lcsim.contrib.JanStrube.tracking.EMap;
-import org.lcsim.contrib.JanStrube.tracking.NewTrack;
-import org.lcsim.recon.vertexing.zvtop4.VectorArithmetic;
-import org.lcsim.spacegeom.SpacePoint;
-
-import Jama.Matrix;
-
-
-/**
- * Class to keep the space-momentum representation of a Track
- * @author jstrube
- *
- */
-public class VertexTrack extends NewTrack {
-
- public VertexTrack(SpacePoint refPoint, EMap parameters, Matrix errorMatrix, int charge, double Bz) {
- super(refPoint, parameters, errorMatrix, charge);
- }
-
- /**
- * Calculates the derivatives of the measurement vector with respect to cartesian spatial coordinates.
- * @param x
- * @param p
- * @param charge
- * @return the derivatives Matrix. Assumes an ordering of ??? in the measurement vector
- */
- public static Matrix getSpatialDerivativeMatrix(SpacePoint x, SpacePoint p, int charge, double Bz) {
- double field = Bz*fieldConversion;
-
- double px0 = p.x() - charge*field;
- double py0 = p.y() + charge*field;
- double pt0 = sqrt(px0*px0 + py0*py0);
-
- // be very explicit about the derivatives
- double dTanLambda_dx = 0;
- double dTanLambda_dy = 0;
- double dTanLambda_dz = 0;
-
- double dPhi_dx = px0*field*charge / (pt0*pt0);
- double dPhi_dy = py0*field*charge / (pt0*pt0);
- double dPhi_dz = 0;
-
- double dOmega_dx = 0;
- double dOmega_dy = 0;
- double dOmega_dz = 0;
-
- double dD0_dx = -py0 / pt0;
- double dD0_dy = px0 / pt0;
- double dD0_dz = 0;
-
- double dZ_dx = -p.z() * px0 / (pt0*pt0);
- double dZ_dy = -p.z() * py0 / (pt0*pt0);
- double dZ_dz = 1;
-
- double[] dh_dx = new double[] {
- dD0_dx, dPhi_dx, dOmega_dx, dZ_dx, dTanLambda_dx
- };
- double[] dh_dy = new double[] {
- dD0_dy, dPhi_dy, dOmega_dy, dZ_dy, dTanLambda_dy
- };
- double[] dh_dz = new double[] {
- dD0_dz, dPhi_dz, dOmega_dz, dZ_dz, dTanLambda_dz
- };
- return new Matrix(new double[][] {dh_dx, dh_dy, dh_dz});
- }
-
-
- public Matrix getMomentumDerivativeMatrix(SpacePoint x, SpacePoint p, int charge, double Bz) {
- double field = Bz*fieldConversion;
- double pt = p.rxy();
- double px0 = p.x() - charge*field;
- double py0 = p.y() + charge*field;
- double pt0 = sqrt(px0*px0 + py0*py0);
- double phi = atan2(p.y(), p.x());
- double phi0 = atan2(py0, px0);
- double l = (phi-phi0)*pt/(charge*field);
-
- double dTanLambda_dpx = -p.z()*p.x() / (pt*pt*pt);
- double dTanLambda_dpy = -p.z()*p.y() / (pt*pt*pt);
- double dTanLambda_dpz = 1/pt;
-
- double dPhi_dpx = -py0 / (pt0*pt0);
- double dPhi_dpy = px0 / (pt0*pt0);
- double dPhi_dpz = 0;
-
- double dOmega_dpx = -charge*field*p.x() / (pt*pt*pt);
- double dOmega_dpy = -charge*field*p.y() / (pt*pt*pt);
- double dOmega_dpz = 0;
-
- double dD0_dpx = -py0/pt0;
- double dD0_dpy = px0/pt0;
- double dD0_dpz = 0;
-
- double dZ_dpx = -p.z() / (field*charge) * (p.y()/(pt*pt) - py0/(pt0*pt0));
- double dZ_dpy = -p.z() / (field*charge) * (p.x()/(pt*pt) - px0/(pt0*pt0));
- double dZ_dpz = -l/pt;
-
- double[] dh_dpx = new double[] {
- dD0_dpx, dPhi_dpx, dOmega_dpx, dZ_dpx, dTanLambda_dpx
- };
- double[] dh_dpy = new double[] {
- dD0_dpy, dPhi_dpy, dOmega_dpy, dZ_dpy, dTanLambda_dpy
- };
- double[] dh_dpz = new double[] {
- dD0_dpz, dPhi_dpz, dOmega_dpz, dZ_dpz, dTanLambda_dpz
- };
- return new Matrix(new double[][] {dh_dpx, dh_dpy, dh_dpz});
- }
-}
lcsim/test/org/lcsim/contrib/JanStrube/vtxFitter
diff -N FitterTest.java
--- FitterTest.java 3 Aug 2006 10:35:54 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,74 +0,0 @@
-package org.lcsim.contrib.JanStrube.vtxFitter;
-
-import static org.lcsim.contrib.JanStrube.vtxFitter.Fitter.add;
-import static org.lcsim.contrib.JanStrube.vtxFitter.Fitter.dot;
-import static org.lcsim.contrib.JanStrube.vtxFitter.Fitter.subtract;
-import junit.framework.TestCase;
-
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.SpacePoint;
-
-import Jama.Matrix;
-
-public class FitterTest extends TestCase {
- Matrix A;
- double[] x;
- double[] y;
- Fitter fitter;
- protected void setUp() throws Exception {
- super.setUp();
- A = new Matrix(2, 3, 2.5);
- x = new double[] {1, 2, 3};
- y = x.clone();
- }
-
- protected void tearDown() throws Exception {
- super.tearDown();
- }
-
-// /*
-// * Test method for 'vtxFitter.Fitter.multiply(Matrix, double[])'
-// */
-// public void testMultiply() {
-// double[] r = multiply(A, x);
-// assertEquals(r.length, A.getRowDimension());
-// assertEquals(r[0], 15.);
-// assertEquals(r[1], 15.);
-// }
-
- /*
- * Test method for 'vtxFitter.Fitter.minus(double[], double[])'
- */
- public void testMinus() {
- double[] r = subtract(x, y);
- assertEquals(r.length, x.length);
- assertEquals(r[0], 0.);
- assertEquals(r[1], 0.);
- assertEquals(r[2], 0.);
- }
-
- /*
- * Test method for 'vtxFitter.Fitter.plus(double[], double[])'
- */
- public void testPlus() {
- double[] r = add(x, y);
- assertEquals(r.length, x.length);
- assertEquals(r[0], 2.);
- assertEquals(r[1], 4.);
- assertEquals(r[2], 6.);
- }
-
- public void testDot() {
- double r = dot(x, y);
- assertEquals(r, 14.);
- }
-
- public void testFit() {
- SpacePoint mom = new CartesianPoint(1, 1, 1);
- SpacePoint pos = new CartesianPoint(0.3, 0.5, 0.2);
- VertexFactory fac = new VertexFactory();
- Vertex particle = fac.makeVertex(pos, mom, 4);
- fitter = new Fitter(particle);
- fitter.fit();
- }
-}
lcsim/test/org/lcsim/contrib/JanStrube/vtxFitter
diff -N VertexFactoryTest.java
--- VertexFactoryTest.java 3 Aug 2006 10:35:54 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,53 +0,0 @@
-package org.lcsim.contrib.JanStrube.vtxFitter;
-
-import org.lcsim.event.ReconstructedParticle;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.SpacePoint;
-
-import junit.framework.TestCase;
-
-public class VertexFactoryTest extends TestCase {
- SpacePoint x;
- SpacePoint p;
- VertexFactory fac;
- protected void setUp() throws Exception {
- super.setUp();
- x = new CartesianPoint(3.6, -4.3, 2.1);
- p = new CartesianPoint(4.2, 2.7, 2.2);
- fac = new VertexFactory();
- }
-
- protected void tearDown() throws Exception {
- super.tearDown();
- }
-
- /*
- * Test method for 'vtxFitter.VertexFactory.makeVertex(SpacePoint, SpacePoint, int)'
- */
- public void testMakeVertex() {
- // positively charged Vertex
- Vertex part = fac.makeVertex(x, p, 1);
- System.out.println(part.origin().toString());
- assertEquals(x.equals(part.origin()), true);
- System.out.println(part.getMomentum().toString());
-// assertEquals(p.equals(part.getMomentum(), 5e-2), true);
-// assertEquals(part.getCharge(), -1.);
-// assertEquals(part.getTracks().size(), 1);
-
- // neutral Vertex
- part = fac.makeVertex(x, p, 4);
-// assertEquals(p.equals(part.getMomentum(), 5e-2), true);
-// assertEquals(x.equals(part.origin()), true);
-// assertEquals(part.getCharge(), 0.);
-// assertEquals(part.getTracks().size(), 4);
-
- // negatively charged Vertex
- part = fac.makeVertex(x, p, 7);
-// assertEquals(p.equals(part.getMomentum(), 5e-2), true);
-// assertEquals(x.equals(part.origin()), true);
-// assertEquals(part.getCharge(), 1.);
-// assertEquals(part.getTracks().size(), 7);
-
- }
-
-}
lcsim/test/org/lcsim/contrib/JanStrube/tracking
diff -N FitterTest.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FitterTest.java 7 Aug 2006 21:30:31 -0000 1.1
@@ -0,0 +1,56 @@
+package org.lcsim.contrib.JanStrube.tracking;
+
+
+import java.util.Random;
+
+import junit.framework.TestCase;
+
+import org.lcsim.contrib.JanStrube.tracking.NewFastMCTrackFactory;
+import org.lcsim.contrib.JanStrube.vtxFitter.Fitter;
+import org.lcsim.contrib.JanStrube.vtxFitter.Vertex;
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.SpacePoint;
+
+import Jama.Matrix;
+
+public class FitterTest extends TestCase {
+ Matrix A;
+ double[] x;
+ double[] y;
+ protected void setUp() throws Exception {
+ super.setUp();
+ A = new Matrix(2, 3, 2.5);
+ x = new double[] {1, 2, 3};
+ y = x.clone();
+ }
+
+ protected void tearDown() throws Exception {
+ super.tearDown();
+ }
+
+// /*
+// * Test method for 'vtxFitter.Fitter.multiply(Matrix, double[])'
+// */
+// public void testMultiply() {
+// double[] r = multiply(A, x);
+// assertEquals(r.length, A.getRowDimension());
+// assertEquals(r[0], 15.);
+// assertEquals(r[1], 15.);
+// }
+
+
+ public void testFit() {
+ Vertex particle = new Vertex();
+ NewFastMCTrackFactory fac = new NewFastMCTrackFactory("sidaug05", 5., false);
+ SpacePoint pos = new CartesianPoint(0, 0, 0);
+ SpacePoint mom1 = new CartesianPoint(1, 2, 3);
+ SpacePoint mom2 = new CartesianPoint(2, -2, 1.5);
+ SpacePoint mom3 = new CartesianPoint(-1, 2, 2);
+ particle.addTrack(fac.getTrack(mom1, pos, pos, -1, new Random(), false), 0);
+ particle.addTrack(fac.getTrack(mom2, pos, pos, 1, new Random(), false), 0);
+ particle.addTrack(fac.getTrack(mom3, pos, pos, -1, new Random(), false), 0);
+ particle.setOrigin(pos);
+ Fitter f = new Fitter(particle, 5.0);
+ f.fit();
+ }
+}
lcsim/test/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.5 -r1.6
--- HelixSwimmerTrackConsistencyTest.java 24 Jul 2006 23:45:41 -0000 1.5
+++ HelixSwimmerTrackConsistencyTest.java 7 Aug 2006 21:30:31 -0000 1.6
@@ -1,5 +1,5 @@
/**
- * @version $Id: HelixSwimmerTrackConsistencyTest.java,v 1.5 2006/07/24 23:45:41 jstrube Exp $
+ * @version $Id: HelixSwimmerTrackConsistencyTest.java,v 1.6 2006/08/07 21:30:31 jstrube Exp $
*/
package org.lcsim.contrib.JanStrube.tracking;
@@ -51,16 +51,16 @@
// swim to pocas of several points and compare position and momentum at these points
double alpha = swimmerRaw.getTrackLengthToPoint(location);
double beta = swimmerTrack.getTrackLengthToPoint(location);
- assertEquals(swimmerRaw.getPointAtDistance(alpha), swimmerTrack.getPointAtDistance(beta), 1e-10);
- assertEquals(swimmerRaw.getMomentumAtDistance(alpha), swimmerTrack.getMomentumAtDistance(beta), 1e-10);
+ assertEquals(swimmerRaw.getPointAtLength(alpha), swimmerTrack.getPointAtLength(beta), 1e-10);
+ assertEquals(swimmerRaw.getMomentumAtLength(alpha), swimmerTrack.getMomentumAtLength(beta), 1e-10);
// use the momentum components as another random point to swim to
alpha = swimmerRaw.getTrackLengthToPoint(momentum);
beta = swimmerTrack.getTrackLengthToPoint(momentum);
// TODO there are obviously some floating point precision issues, probably related to the trigonometric functions.
// Investigation can be put off to a later date
- assertEquals(swimmerRaw.getPointAtDistance(alpha), swimmerTrack.getPointAtDistance(beta), 1e-10);
- assertEquals(swimmerRaw.getMomentumAtDistance(alpha), swimmerTrack.getMomentumAtDistance(beta), 1e-10);
+ assertEquals(swimmerRaw.getPointAtLength(alpha), swimmerTrack.getPointAtLength(beta), 1e-10);
+ assertEquals(swimmerRaw.getMomentumAtLength(alpha), swimmerTrack.getMomentumAtLength(beta), 1e-10);
}
public void testHelix_TrackFactoryConsistency2() {
@@ -88,18 +88,18 @@
double alpha_track = swimmerTrack.getTrackLengthToPoint(p1)+alpha;
assertEquals(alpha_raw, alpha_track, 1e-10);
// then the POCA to that point must be the same, too
- assertEquals(swimmerRaw.getPointAtDistance(alpha_raw), swimmerTrack.getPointAtDistance(alpha_track-alpha), 1e-10);
+ assertEquals(swimmerRaw.getPointAtLength(alpha_raw), swimmerTrack.getPointAtLength(alpha_track-alpha), 1e-10);
// d0 is the distance between P0 and Pref in the x-y plane
double d0 = t.getParameter(Track.ParameterName.d0);
- SpacePoint origin = swimmerTrack.getPointAtDistance(0);
+ SpacePoint origin = swimmerTrack.getPointAtLength(0);
double dist = VectorArithmetic.subtract(origin, t.getReferencePoint()).rxy();
assertEquals(abs(d0), abs(dist), 1e-6);
// The points at these two locations should be the same, too.
for (int i=-70; i<70; i++) {
double alpha_i = i/10;
- assertEquals(swimmerRaw.getPointAtDistance(alpha_i), swimmerTrack.getPointAtDistance(alpha_i-alpha), 1e-10);
+ assertEquals(swimmerRaw.getPointAtLength(alpha_i), swimmerTrack.getPointAtLength(alpha_i-alpha), 1e-10);
}
}
@@ -130,8 +130,8 @@
NewTrack t1 = new NewTrack(referencePoint, tp, new Matrix(0, 0), charge);
swimmerRaw.setTrack(t);
swimmerTrack.setTrack(t1);
- assertEquals(swimmerRaw.getMomentumAtDistance(0), swimmerTrack.getMomentumAtDistance(0));
- assertEquals(swimmerRaw.getPointAtDistance(0), swimmerTrack.getPointAtDistance(0));
+ assertEquals(swimmerRaw.getMomentumAtLength(0), swimmerTrack.getMomentumAtLength(0));
+ assertEquals(swimmerRaw.getPointAtLength(0), swimmerTrack.getPointAtLength(0));
}
/* (non-Javadoc)
lcsim/src/org/lcsim/recon/vertexing/zvtop4
diff -u -r1.12 -r1.13
--- VectorArithmetic.java 18 Jul 2006 03:26:54 -0000 1.12
+++ VectorArithmetic.java 7 Aug 2006 21:30:32 -0000 1.13
@@ -12,7 +12,7 @@
/**
* Class to perform basic vector arithmetic on Hep3Vectors
* @author jstrube
- * @version $Id: VectorArithmetic.java,v 1.12 2006/07/18 03:26:54 jstrube Exp $
+ * @version $Id: VectorArithmetic.java,v 1.13 2006/08/07 21:30:32 jstrube Exp $
*/
final public class VectorArithmetic {
@@ -32,6 +32,15 @@
return a.x() * b.x() + a.y() * b.y() + a.z() * b.z();
}
+ public static double dot(double[] a, double[] b) {
+ if (a.length != b.length)
+ throw new IllegalArgumentException("dimensions don't match");
+ double result = 0;
+ for (int i=0; i<a.length; i++)
+ result += a[i]*b[i];
+ return result;
+ }
+
public static Hep3Vector cross(Hep3Vector a, Hep3Vector b) {
double x = a.y() * b.z() - a.z() * b.y();
double y = a.z() * b.x() - a.y() * b.x();
@@ -80,7 +89,13 @@
}
public static double[] subtract(double[] a, double[] b) {
- return new double[] { a[0] - b[0], a[1] - b[1], a[2] - b[2] };
+ if (a.length != b.length)
+ throw new IllegalArgumentException("dimensions do not match");
+ double[] result = new double[a.length];
+ for (int i=0; i<a.length; i++) {
+ result[i] = a[i] - b[i];
+ }
+ return result;
};
/**
@@ -126,7 +141,16 @@
return new BasicHep3Vector(x);
}
-
+ public static double[] add(double[] a, double[] b) {
+ if (a.length != b.length)
+ throw new IllegalArgumentException("dimensions do not match");
+ double[] result = new double[a.length];
+ for (int i=0; i<a.length; i++) {
+ result[i] = a[i] + b[i];
+ }
+ return result;
+ }
+
/**
* Returns the distance between two space points.
* @param spt1 SpacePoint 1
CVSspam 0.2.8