2 added + 6 modified, total 8 files
lcsim/src/org/lcsim/util/swim
diff -u -r1.20 -r1.21
--- Helix.java 30 May 2006 06:17:50 -0000 1.20
+++ Helix.java 4 Jul 2006 01:41:47 -0000 1.21
@@ -21,7 +21,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.20 2006/05/30 06:17:50 jstrube Exp $
+ * @version $Id: Helix.java,v 1.21 2006/07/04 01:41:47 jstrube Exp $
*/
public class Helix implements Trajectory
{
@@ -51,7 +51,7 @@
sinLambda = sin(lambda);
xCenter = origin.x() + radius*sin(phi);
yCenter = origin.y() - radius*cos(phi);
- setSpacialParameters();
+ setSpatialParameters();
}
/**
@@ -182,7 +182,7 @@
* Sets the parametereization in terms of "momentum" and charge
*
*/
- private void setSpacialParameters() {
+ private void setSpatialParameters() {
abs_r = abs(radius);
px = abs_r*cosPhi;
py = abs_r*sinPhi;
lcsim/src/org/lcsim/util/swim
diff -u -r1.14 -r1.15
--- HelixSwimmer.java 28 Jun 2006 01:54:50 -0000 1.14
+++ HelixSwimmer.java 4 Jul 2006 01:41:47 -0000 1.15
@@ -14,17 +14,17 @@
* For more info on swimming see <a href="doc-files/transport.pdf">this paper</a>
* by Paul Avery.
* @author tonyj
- * @version $Id: HelixSwimmer.java,v 1.14 2006/06/28 01:54:50 jstrube Exp $
+ * @version $Id: HelixSwimmer.java,v 1.15 2006/07/04 01:41:47 jstrube Exp $
*/
public class HelixSwimmer
{
- class SpatialParameters
+ public class SpatialParameters
{
- double px;
- double py;
- double pz;
- int charge;
- boolean isInvalid = true;
+ public double px;
+ public double py;
+ public double pz;
+ public int charge;
+ public boolean isInvalid = true;
}
private double field;
private Trajectory trajectory;
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -N VertexTrackTest.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ VertexTrackTest.java 4 Jul 2006 01:41:48 -0000 1.1
@@ -0,0 +1,52 @@
+package org.lcsim.contrib.JanStrube.vtxFitter;
+
+import hep.physics.vec.Hep3Vector;
+
+import org.lcsim.event.Track;
+import org.lcsim.mc.fast.tracking.MCFastTrackFactory;
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.util.swim.Helix;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.util.swim.HelixSwimmer.SpatialParameters;
+
+import junit.framework.TestCase;
+
+public class VertexTrackTest extends TestCase {
+
+ HelixSwimmer helix;
+ MCFastTrackFactory tf;
+ VertexTrack vtxTrack;
+ protected void setUp() throws Exception {
+ super.setUp();
+ helix = new HelixSwimmer(5);
+ tf = new MCFastTrackFactory();
+ SpacePoint momentum = new CartesianPoint(1, 2, 3);
+ SpacePoint location = new CartesianPoint(1, 2, 3);
+ int charge = -1;
+ Track t = tf.getMCTrack(momentum, location, charge);
+ helix.setTrack(t);
+ vtxTrack = new VertexTrack(t);
+
+ }
+
+ protected void tearDown() throws Exception {
+ super.tearDown();
+ }
+
+ public void testGetPointAtLength() {
+ Hep3Vector p1 = helix.getPointAtDistance(0.5);
+ double[] p2 = vtxTrack.getPointAtLength(0.5);
+ assertEquals(p1.x(), p2[0], 1e-10);
+ assertEquals(p1.y(), p2[1], 1e-10);
+ assertEquals(p1.z(), p2[2], 1e-10);
+ }
+
+ public void testGetMomentumAtLength() {
+ SpatialParameters mom1 = helix.getSpatialParameters();
+ double[] mom2 = vtxTrack.getMomentumAtLength(0);
+ assertEquals(mom1.px, mom2[0], 1e-10);
+ assertEquals(mom1.py, mom2[1], 1e-10);
+ assertEquals(mom1.pz, mom2[2], 1e-10);
+ }
+}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -N VertexTrack.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ VertexTrack.java 4 Jul 2006 01:41:48 -0000 1.1
@@ -0,0 +1,302 @@
+/**
+ * @version $Id: VertexTrack.java,v 1.1 2006/07/04 01:41:48 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 org.lcsim.constants.Constants.fieldConversion;
+
+import java.util.List;
+
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+
+import Jama.Matrix;
+
+
+/**
+ * Class to keep the space-momentum representation of a Track
+ * @author jstrube
+ *
+ */
+public class VertexTrack implements Track {
+ private Track _track;
+ private static double Bz = 5.;
+ /**
+ * Creates a new instance from an existing Track
+ */
+ public VertexTrack(Track t) {
+ _track = t;
+ }
+
+ public double[] getPointAtLength(double s) {
+ double phi0 = _track.getTrackParameter(Track.ParameterName.phi0.ordinal());
+ double omega = _track.getTrackParameter(Track.ParameterName.omega.ordinal());
+ double r = 1/omega;
+ double d0 = _track.getTrackParameter(Track.ParameterName.d0.ordinal());
+ double z0 = _track.getTrackParameter(Track.ParameterName.z0.ordinal());
+ double tanLambda = _track.getTrackParameter(Track.ParameterName.s.ordinal());
+ double phi = phi0 + omega*s;
+
+ return new double[] {
+ r*sin(phi) - (r+d0)*sin(phi0)
+ , -r*cos(phi) + (r+d0)*cos(phi0)
+ , z0 + s*tanLambda
+ };
+ }
+
+ public double[] getMomentumAtLength(double s) {
+ double phi0 = _track.getTrackParameter(Track.ParameterName.phi0.ordinal());
+ double omega = _track.getTrackParameter(Track.ParameterName.omega.ordinal());
+ double r = 1/omega;
+ double tanLambda = _track.getTrackParameter(Track.ParameterName.s.ordinal());
+ int charge = _track.getCharge();
+ double phi = phi0 + omega*s;
+ double pt = charge*fieldConversion*r;
+ return new double[] {
+ pt*cos(phi), pt*sin(phi), pt*tanLambda
+ };
+ }
+
+ /**
+ * Calculates the 5-parameter measurement vector from position and momentum vectors
+ * @param x Cartesian position vector (3 parameters)
+ * @param p Cartesian momentum vector (3 parameters)
+ * @return the 5-parameter vector with the parameters as specified by {@link ParameterName#values()}
+ */
+ public static double[] getMeasurementVector(double[] x, double[] p, int charge) {
+ double pt = sqrt(p[0]*p[0] + p[1]*p[1]);
+ double px0 = p[0] + charge*Bz*fieldConversion;
+ double py0 = p[1] - charge*Bz*fieldConversion;
+ double pt0 = sqrt(px0*px0 + py0*py0);
+ double phi = atan2(p[1], p[0]);
+ double phi0 = atan2(py0, px0);
+ double l = (phi-phi0)*pt/(charge*fieldConversion);
+
+ double d0 = (pt0-pt)/(charge*fieldConversion);
+ // FIXME signed/unsigned ? BaBar convention vs. org.lcsim convention
+ double omega = charge*fieldConversion/pt;
+ double z0 = x[2] - l*p[2]/pt;
+ double tanLambda = p[2]/pt;
+ return new double[]{d0, phi0, omega, tanLambda, z0};
+ }
+
+ public static Matrix getSpatialDerivativeMatrix(double[] x, double[] p, int charge) {
+ double px0 = p[0] + charge*Bz*fieldConversion;
+ double py0 = p[1] - charge*Bz*fieldConversion;
+ double pt0 = sqrt(px0*px0 + py0*py0);
+
+ double[] dh_dx = new double[] {
+ -py0/pt0, -charge*fieldConversion*px0/(pt0*pt0), 0, -p[2]*px0/(pt0*pt0), 0
+ };
+ double[] dh_dy = new double[] {
+ px0/pt0, -charge*fieldConversion*py0/(pt0*pt0), 0, -p[2]*px0/(pt0*pt0), 0
+ };
+ double[] dh_dz = new double[] {
+ 0, 0, 0, 1, 0
+ };
+ return new Matrix(new double[][] {dh_dx, dh_dy, dh_dz});
+ }
+
+
+ public Matrix getMomentumDerivativeMatrix(double[] x, double[] p, int charge) {
+ double pt = sqrt(p[0]*p[0] + p[1]*p[1]);
+ double px0 = p[0] + charge*Bz*fieldConversion;
+ double py0 = p[1] - charge*Bz*fieldConversion;
+ double pt0 = sqrt(px0*px0 + py0*py0);
+ double phi = atan2(p[1], p[0]);
+ double phi0 = atan2(py0, px0);
+ double l = (phi-phi0)*pt/(charge*fieldConversion);
+
+ double[] dh_dpx = new double[] {
+ 1/(charge*fieldConversion)*(px0/pt0 - p[0]/pt)
+ , -py0/(pt0*pt0)
+ , -charge*fieldConversion*p[0]/(pt*pt*pt)
+ , -p[2]/(charge*fieldConversion)*(py0/(pt0*pt0)-p[1]/(pt0*pt0))
+ , -p[2]*p[0]/(pt*pt*pt)
+ };
+ double[] dh_dpy = new double[] {
+ 1/(charge*fieldConversion)*(py0/pt0 - p[1]/pt)
+ , -px0/(pt0*pt0)
+ , -charge*fieldConversion*p[1]/(pt*pt*pt)
+ , -p[2]/(charge*fieldConversion)*(px0/(pt0*pt0)-p[0]/(pt0*pt0))
+ , -p[2]*p[1]/(pt*pt*pt)
+ };
+ double[] dh_dpz = new double[] {
+ 0, 0, 0, -l/pt, 1/pt
+ };
+ return new Matrix(new double[][] {dh_dpx, dh_dpy, dh_dpz});
+ }
+
+ /**
+ * @see org.lcsim.event.Track#fitSuccess()
+ */
+ public boolean fitSuccess() {
+ return _track.fitSuccess();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getCharge()
+ */
+ public int getCharge() {
+ return _track.getCharge();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getChi2()
+ */
+ public double getChi2() {
+ return _track.getChi2();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getErrorMatrix()
+ */
+ public double[][] getErrorMatrix() {
+ return _track.getErrorMatrix();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getErrorMatrixElement(int, int)
+ */
+ public double getErrorMatrixElement(int i, int j) {
+ return _track.getErrorMatrixElement(i, j);
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getMomentum()
+ */
+ public double[] getMomentum() {
+ return _track.getMomentum();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getNDF()
+ */
+ public int getNDF() {
+ return _track.getNDF();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getPX()
+ */
+ public double getPX() {
+ return _track.getPX();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getPY()
+ */
+ public double getPY() {
+ return _track.getPY();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getPZ()
+ */
+ public double getPZ() {
+ return _track.getPZ();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getRadiusOfInnermostHit()
+ */
+ public double getRadiusOfInnermostHit() {
+ return _track.getRadiusOfInnermostHit();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getReferencePoint()
+ */
+ public double[] getReferencePoint() {
+ return _track.getReferencePoint();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getReferencePointX()
+ */
+ public double getReferencePointX() {
+ return _track.getReferencePointX();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getReferencePointY()
+ */
+ public double getReferencePointY() {
+ return _track.getReferencePointY();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getReferencePointZ()
+ */
+ public double getReferencePointZ() {
+ return _track.getReferencePointZ();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getSubdetectorHitNumbers()
+ */
+ public int[] getSubdetectorHitNumbers() {
+ return _track.getSubdetectorHitNumbers();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getTrackParameter(int)
+ */
+ public double getTrackParameter(int i) {
+ return _track.getTrackParameter(i);
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getTrackParameters()
+ */
+ public double[] getTrackParameters() {
+ return _track.getTrackParameters();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getTrackerHits()
+ */
+ public List<TrackerHit> getTrackerHits() {
+ return _track.getTrackerHits();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getTracks()
+ */
+ public List<Track> getTracks() {
+ return _track.getTracks();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getType()
+ */
+ public int getType() {
+ return _track.getType();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getdEdx()
+ */
+ public double getdEdx() {
+ return _track.getdEdx();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#getdEdxError()
+ */
+ public double getdEdxError() {
+ return _track.getdEdxError();
+ }
+
+ /**
+ * @see org.lcsim.event.Track#isReferencePointPCA()
+ */
+ public boolean isReferencePointPCA() {
+ return _track.isReferencePointPCA();
+ }
+
+}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -u -r1.3 -r1.4
--- VertexFactory.java 28 Jun 2006 20:01:18 -0000 1.3
+++ VertexFactory.java 4 Jul 2006 01:41:48 -0000 1.4
@@ -13,10 +13,13 @@
/**
* Class to create ReconstructedParticle objects consisting
- * of tracks with random (within limits) properties
+ * 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;
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -u -r1.4 -r1.5
--- Fitter.java 28 Jun 2006 04:48:33 -0000 1.4
+++ Fitter.java 4 Jul 2006 01:41:48 -0000 1.5
@@ -1,6 +1,6 @@
package org.lcsim.contrib.JanStrube.vtxFitter;
/**
- * @version $Id: Fitter.java,v 1.4 2006/06/28 04:48:33 jstrube Exp $
+ * @version $Id: Fitter.java,v 1.5 2006/07/04 01:41:48 jstrube Exp $
*/
import org.lcsim.event.ReconstructedParticle;
@@ -9,7 +9,6 @@
import Jama.Matrix;
public class Fitter {
-
ReconstructedParticle particle;
// convergence criterion
double deltaChi2 = 1e-4;
@@ -28,20 +27,18 @@
Matrix V_prev;
// inverse of covariance matrix
Matrix G_prev;
- Matrix G_B;
// linearity factors of the measeurement eq.
// in x
Matrix A_k, A_prev, A;
// and q
Matrix B_k, B_prev, B;
- Matrix W_k;
// residual
double[] c_k0;
// chi2
- double chi2, chi2_prev;
+ double chi2_prev;
double[] r_k;
@@ -81,29 +78,33 @@
smoothe(t);
}
}
- private void filter(Track t) {
+ 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))));
+
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();
- W_k = (B_k.transpose().times(G_k.times(B_k))).times(-1);
- G_B = G_k.minus((G_k.times(B_k.times(W_k.times(B_k.transpose().times(G_k))))));
- 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 + dot(r_k, G_k.times(r_k)) + dot(subtract(x_k, x_prev), C_prev.inverse().times(subtract(x_k, x_prev)));
r_k = subtract(m_k, add(c_k0, add(A_k.times(x_k), B_k.times(q_k))));
C_prev = C_k;
+ 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))));
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)))))))));
@@ -153,6 +154,16 @@
static Matrix add(Matrix a, Matrix b) {
return a.plus(b);
}
+
+ // matrix A_c in the Grab and Luchsinger paper
+ static Matrix chargedSpatialLinearity(Track t) {
+
+ return new Matrix(0, 0);
+ }
+
+ static Matrix chargedMomentumLinearity(Track t) {
+ return new Matrix(0, 0);
+ }
static double dot(double[] a, double[] b) {
if (a.length != b.length)
@@ -163,4 +174,11 @@
return result;
}
+ static Matrix calcSpatialLinearityMatrix(Track t) {
+ return new Matrix(5, 5);
+ }
+
+ static Matrix calcMomentumLinearityMatrix(Track t) {
+ return new Matrix(5, 5);
+ }
}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -u -r1.3 -r1.4
--- FitterTest.java 27 Jun 2006 09:00:14 -0000 1.3
+++ FitterTest.java 4 Jul 2006 01:41:48 -0000 1.4
@@ -18,14 +18,12 @@
Matrix A;
double[] x;
double[] y;
- ReconstructedParticle part;
+ 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();
- part = new BaseReconstructedParticle();
- Fitter fitter = new Fitter(part);
}
protected void tearDown() throws Exception {
@@ -72,9 +70,9 @@
public void testFit() {
SpacePoint mom = new CartesianPoint(1, 1, 1);
SpacePoint pos = new CartesianPoint(0.3, 0.5, 0.2);
- MCFastTrackFactory fac = new MCFastTrackFactory();
- Track t = fac.getMCTrack(mom, pos, -1);
- part.addTrack(t);
-
+ VertexFactory fac = new VertexFactory();
+ Vertex particle = fac.makeVertex(pos, mom, 4);
+ fitter = new Fitter(particle);
+ fitter.fit();
}
}
lcsim/src/org/lcsim/event/base
diff -u -r1.5 -r1.6
--- BaseTrack.java 29 Mar 2006 18:51:26 -0000 1.5
+++ BaseTrack.java 4 Jul 2006 01:41:49 -0000 1.6
@@ -3,7 +3,7 @@
*
* Created on March 24, 2006, 9:18 PM
*
- * $Id: BaseTrack.java,v 1.5 2006/03/29 18:51:26 ngraf Exp $
+ * $Id: BaseTrack.java,v 1.6 2006/07/04 01:41:49 jstrube Exp $
*/
package org.lcsim.event.base;
@@ -38,6 +38,12 @@
protected List<TrackerHit> _hits;
protected int _type;
+ protected static final int D0 = Track.ParameterName.d0.ordinal();
+ protected static final int PHI = Track.ParameterName.phi0.ordinal();
+ protected static final int OMEGA = Track.ParameterName.omega.ordinal();
+ protected static final int TANLAMBDA = Track.ParameterName.s.ordinal();
+ protected static final int Z0 = Track.ParameterName.z0.ordinal();
+
/** Creates a new instance of BaseTrack */
public BaseTrack()
{
@@ -84,16 +90,15 @@
* @see Track.Parameters
* @param params The array of track parameters.
*/
- // TODO remove hard-coded track parameter indices.
public void setTrackParameters(double[] params, double magneticField)
{
System.arraycopy(params, 0, _parameters, 0, 5);
- double omega = _parameters[2];
+ double omega = _parameters[OMEGA];
if(abs(omega) < 0.0000001) omega=0.0000001;
double Pt = abs((1./omega) * magneticField* Constants.fieldConversion);
- _momentum[0] = Pt * Math.cos(_parameters[1]);
- _momentum[1] = Pt * Math.sin(_parameters[1]);
- _momentum[2] = Pt * _parameters[4];
+ _momentum[0] = Pt * Math.cos(_parameters[PHI]);
+ _momentum[1] = Pt * Math.sin(_parameters[PHI]);
+ _momentum[2] = Pt * _parameters[TANLAMBDA];
_charge = (int) signum(omega);
}
@@ -179,11 +184,11 @@
int lastDot = className.lastIndexOf('.');
if(lastDot!=-1)className = className.substring(lastDot+1);
StringBuffer sb = new StringBuffer(className+": Type: "+_type+" charge: "+ _charge+"\n");
- sb.append("d0= "+_parameters[0]+"\n");
- sb.append("phi0= "+_parameters[1]+"\n");
- sb.append("curvature: "+_parameters[2]+"\n");
- sb.append("z0= "+_parameters[3]+"\n");
- sb.append("tanLambda= "+_parameters[4]+"\n");
+ sb.append("d0= "+_parameters[D0]+"\n");
+ sb.append("phi0= "+_parameters[PHI]+"\n");
+ sb.append("curvature: "+_parameters[OMEGA]+"\n");
+ sb.append("z0= "+_parameters[Z0]+"\n");
+ sb.append("tanLambda= "+_parameters[TANLAMBDA]+"\n");
sb.append(" px="+getPX()+" py= "+getPY()+" pz= "+getPZ());
return sb.toString();
}
CVSspam 0.2.8
