Print

Print
Commit in lcsim/src/org/lcsim/recon/vertexing/zvtop4 on MAIN
ZvUtil.java+25-121.7 -> 1.8
ZvTrack.java+270-511.14 -> 1.15
ZvTop.java+125-121.9 -> 1.10
ZvMaximum.java+8-31.1 -> 1.2
+428-78
4 modified files
- changed transformLocation() to rotateLocation()
- stupidTwoTrackMax was a little too stupid, instead of looking at random points in space, it now looks at 100 points between one of 1000 points on one track and the POCA on the other track. Formula derived with help of Maple. Still needs debugging.

In principle the formula for the closes distance between two line should be deterministic, no ?
lcsim/src/org/lcsim/recon/vertexing/zvtop4
ZvUtil.java 1.7 -> 1.8
diff -u -r1.7 -r1.8
--- ZvUtil.java	30 Jul 2005 19:22:26 -0000	1.7
+++ ZvUtil.java	1 Aug 2005 18:09:50 -0000	1.8
@@ -21,6 +21,8 @@

 import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.multiply;
 import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.unit;
 

+import static java.lang.Math.sqrt;
+

 /**
  * Repository of utility functions
  * @author jstrube

@@ -49,25 +51,36 @@

     private static SpacePoint stupidTwoTrackMax(ZvTrack iTrack, ZvTrack jTrack) {
         SpacePoint maxLocation = new SpacePoint();
         // limit the region where to look to the vertex detector.

-        double maxX = 100;
-        double maxY = 100;
-        double maxZ = 500;

+        // FIXME look at the real dimensions
+        final double maxX = 100;
+        final double maxY = 100;
+        final double maxZ = 500;
+        final int maxStep = 1000;

         double maxOverlap = 0;

-        for (int x=1; x<=100; ++x) {
-            for (int y=1; y<=100; ++y) {
-                for (int z=1; z<=50; ++z) {
-                    SpacePoint location = new CartesianPoint(maxX/x, maxY/y, maxZ/z);
-                    double overlap = iTrack.getTubeValue(location)*jTrack.getTubeValue(location);
-                    if (overlap > maxOverlap) {
-                        maxOverlap = overlap;
-                        maxLocation = location;
-                    }

+        // cut the Track into maxStep steps and calculate the distance
+        // to the other parabola
+        for (int s=0; s <= maxStep; ++s) {
+            double lengthParameter = - sqrt(maxX*maxX + maxY*maxY) + s*maxY/maxStep;
+            SpacePoint thisITrackPoint = iTrack.getTubeCenter(lengthParameter);
+            double closestJTrackLength = jTrack.closestParabolaPoint(thisITrackPoint);
+            SpacePoint closestJTrackPoint = jTrack.getTubeCenter(closestJTrackLength);
+            SpacePath t = new SpacePath();
+            t.setOrigin(thisITrackPoint);
+            t.setEndPoint(closestJTrackPoint);
+            for(int iPoint=1; iPoint<=100; ++iPoint) {
+                SpacePoint checkPoint = t.getPointAtLength(1/iPoint);
+                double overlap = iTrack.getTubeValue(checkPoint) * jTrack.getTubeValue(checkPoint); 
+                if (overlap > maxOverlap) {
+                    maxOverlap = overlap;
+                    maxLocation = checkPoint;

                 }
             }
         }
         return maxLocation;
     }
     

+    
+    

     /**
      * Sums the momenta of all tracks in the list
      * @param list
lcsim/src/org/lcsim/recon/vertexing/zvtop4
ZvTrack.java 1.14 -> 1.15
diff -u -r1.14 -r1.15
--- ZvTrack.java	30 Jul 2005 19:22:26 -0000	1.14
+++ ZvTrack.java	1 Aug 2005 18:09:50 -0000	1.15
@@ -10,6 +10,7 @@

 
 import org.lcsim.event.Track;
 import org.lcsim.event.TrackerHit;

+import org.lcsim.geometry.compact.converter.lcdd.util.Rotation;

 import org.lcsim.spacegeom.CartesianPoint;
 import org.lcsim.spacegeom.SpacePoint;
 

@@ -21,12 +22,14 @@

 
 /**
  * Representation of the Gaussian tubes

+ * 

  * @author jstrube
  */
 public class ZvTrack implements Track {
 

-    // TODO ksi and kappa need to be calculated -- keep in mind the charge of the track
-    

+    // TODO ksi and kappa need to be calculated -- keep in mind the charge of
+    // the track
+

     SpacePoint _referencePoint;
     Hep3Vector _momentum;
     ZvParameters _parameters;

@@ -40,6 +43,7 @@

     Track orgTrack;
 
     double trackRotationAngle;

+

     // package wide constructor ensures I don't have to worry about anybody
     // creating invalid ZvTracks.
     // Only used to create Track-IP

@@ -52,7 +56,7 @@

         fitSuccess = false;
         updateMomentum();
     }

-    

+

     // only useful for package testing
     // NEVER use this as a user, as nothing guarantees consistency
     ZvTrack(ZvParameters params, SpacePoint refPoint, double[][] errorMatrix) {

@@ -64,6 +68,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getCharge()
      */
     public int getCharge() {

@@ -72,14 +77,17 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getReferencePoint()
      */

-    @Deprecated public double[] getReferencePoint() {

+    @Deprecated
+    public double[] getReferencePoint() {

         return _referencePoint.getCartesianArray();
     }
 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getReferencePointX()
      */
     public double getReferencePointX() {

@@ -88,6 +96,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getReferencePointY()
      */
     public double getReferencePointY() {

@@ -96,6 +105,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getReferencePointZ()
      */
     public double getReferencePointZ() {

@@ -104,6 +114,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#isReferencePointPCA()
      */
     public boolean isReferencePointPCA() {

@@ -112,6 +123,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getMomentum()
      */
     public double[] getMomentum() {

@@ -120,6 +132,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getPX()
      */
     public double getPX() {

@@ -128,6 +141,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getPY()
      */
     public double getPY() {

@@ -136,6 +150,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getPZ()
      */
     public double getPZ() {

@@ -144,15 +159,16 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#fitSuccess()
      */
     public boolean fitSuccess() {
         return fitSuccess;
     }
 

-

     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getTrackParameter(int)
      */
     public double zvParameter(ZvParameterNames name) {

@@ -161,8 +177,10 @@

 
     /**
      * returns the parameter by index

+     * 

      * @deprecated, because the order is arbitrary

-     * @param i <em>Index</em> of the parameter

+     * @param i
+     *            <em>Index</em> of the parameter

      * @return the parameter
      */
     @Deprecated

@@ -170,9 +188,9 @@

         return orgTrack.getTrackParameter(i);
     }
 

-    

     /**
      * returns an array of parameters

+     * 

      * @deprecated, because the order is arbitrary
      * @return the parameters as double[]
      */

@@ -183,6 +201,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getTrackParameters()
      */
     public ZvParameters zvParameters() {

@@ -191,6 +210,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getErrorMatrixElement(int, int)
      */
     public double getZvErrorMatrixElement(int i, int j) {

@@ -199,6 +219,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getErrorMatrix()
      */
     public double[][] getZvErrorMatrix() {

@@ -207,27 +228,33 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getErrorMatrixElement(int, int)
      */
     /**
      * @deprecated because it's expensive
      */

-    @Deprecated public double getErrorMatrixElement(int i, int j) {

+    @Deprecated
+    public double getErrorMatrixElement(int i, int j) {

         return convert2TrackErrorMatrix()[i][j];
     }
 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getErrorMatrix()
      */
     /**
      * @deprecated because it's expensive
      */

-    @Deprecated public double[][] getErrorMatrix() {

+    @Deprecated
+    public double[][] getErrorMatrix() {

         return convert2TrackErrorMatrix();
     }

+

     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getChi2()
      */
     public double getChi2() {

@@ -236,6 +263,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getNDF()
      */
     public int getNDF() {

@@ -244,6 +272,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getdEdx()
      */
     public double getdEdx() {

@@ -252,14 +281,16 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getdEdxError()
      */
     public double getdEdxError() {

-        return orgTrack.getdEdxError(); 

+        return orgTrack.getdEdxError();

     }
 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getRadiusOfInnermostHit()
      */
     public double getRadiusOfInnermostHit() {

@@ -268,6 +299,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getSubdetectorHitNumbers()
      */
     public int[] getSubdetectorHitNumbers() {

@@ -276,6 +308,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getTracks()
      */
     public List<Track> getTracks() {

@@ -284,6 +317,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getTrackerHits()
      */
     public List<TrackerHit> getTrackerHits() {

@@ -292,6 +326,7 @@

 
     /*
      * (non-Javadoc)

+     * 

      * @see org.lcsim.event.Track#getType()
      */
     public int getType() {

@@ -299,72 +334,95 @@

     }
 
     /**

-     * Calculates a Gaussian uncertainty around the track.
-     * A parabolic approximation is made instead of a helical parametrization
-     * The apex is at (x0, 0, z0), which is the POCA in a co-ordinate system that
-     * has been rotated in the x-y plane, so that the track momentum at the POCA
-     * in the x-y plane is parallel to the y-axis. 

+     * Calculates a Gaussian uncertainty around the track. A parabolic
+     * approximation is made instead of a helical parametrization The apex is at
+     * (x0, 0, z0), which is the POCA in a co-ordinate system that has been
+     * rotated in the x-y plane, so that the track momentum at the POCA in the
+     * x-y plane is parallel to the y-axis.

      * 

-     * @param location The spatial position at which to evaluate the "track probability"
-     * in the lab frame

+     * @param location
+     *            The spatial position at which to evaluate the "track
+     *            probability" in the lab frame

      * @return the scalar value at that location
      */
     public double getTubeValue(SpacePoint location) {

-        SpacePoint r = transformLocation(location);
-        SpacePoint refPrime = transformLocation(_referencePoint);
-        double termA = (r.x() - refPrime.x() - _parameters.get(ZvParameterNames.kappa)
-                * r.y() * r.y())
-                / sigmaT;
-        double termB = (r.z() - refPrime.z() - _parameters.get(ZvParameterNames.s)

+        SpacePoint r = rotateLocation(location);
+        SpacePoint refPrime = rotateLocation(_referencePoint);
+        double termA = (r.x() - refPrime.x() - _parameters.get(ZvParameterNames.kappa) * r.y()

                 * r.y())

+                / sigmaT;
+        double termB = (r.z() - refPrime.z() - _parameters.get(ZvParameterNames.s) * r.y())

                 / sigmaL;
         double result = Math.exp(-0.5 * (termA * termA + termB * termB));

-        // System.out.printf("termA: %f\ttermB: %f\tresult: %f", termA, termB, result);

+        // System.out.printf("termA: %f\ttermB: %f\tresult: %f", termA, termB,
+        // result);

         return result;
     }
 
     /**
      * Projection of the parabolic approximation of the track into the x-y plane
      * Useful for visualization

+     * 

      * @see getTubeValue

-     * @param location The spatial position at which to evaluate the "track probability"
-     * in the lab frame

+     * @param location
+     *            The spatial position at which to evaluate the "track
+     *            probability" in the lab frame

      * @return the scalar value at that location
      */
     public double getTubeProjection(SpacePoint location) {

-        SpacePoint r = transformLocation(location);
-        SpacePoint refPrime = transformLocation(_referencePoint);
-        double termA = (r.x() - refPrime.x() - _parameters.get(ZvParameterNames.kappa)
-                * r.y() * r.y())

+        SpacePoint r = rotateLocation(location);
+        SpacePoint refPrime = rotateLocation(_referencePoint);
+        double termA = (r.x() - refPrime.x() - _parameters.get(ZvParameterNames.kappa) * r.y()
+                * r.y())

                 / sigmaT;
         double result = Math.exp(-0.5 * termA * termA);

-        // System.out.printf("termA: %f\ttermB: %f\tresult: %f", termA, termB, result);

+        // System.out.printf("termA: %f\ttermB: %f\tresult: %f", termA, termB,
+        // result);

         return result;
     }
 

-    

+    /**
+     * Returns a point in space in the center of the tube. The parabola is
+     * parametrized by a length parameter s The value for s=0 ist the reference
+     * point
+     * 
+     * @param lengthParameter
+     *            The length parameter of the parabola
+     * @return The point in space at the length lengthParameter
+     */
+    SpacePoint getTubeCenter(double lengthParameter) {
+        SpacePoint startingPoint = rotateLocation(referencePoint());
+        double theta = trackRotationAngle;
+        double yPrime = lengthParameter;
+        double xPrime = startingPoint.x() + _parameters.kappa * yPrime * yPrime;
+        double zPrime = startingPoint.z() + _parameters.s * yPrime;
+        double x = xPrime * cos(theta) - yPrime * sin(theta);
+        double y = yPrime * cos(theta) + xPrime * sin(theta);
+        return new CartesianPoint(x, y, zPrime);
+    }
+

     // TODO implementation
     public ZvFitStatus getSwimStatus() {
         return ZvFitStatus.FAILED_ON_ENTRY;
     }
 

-    // Rotate the coordinate system, so that the y axis is parallel to the momentum vector
-    SpacePoint transformLocation(SpacePoint oldPoint) {

+    // Rotate the coordinate system, so that the y axis is parallel to the
+    // momentum vector
+    SpacePoint rotateLocation(SpacePoint oldPoint) {

         // z stays the same
         double z = oldPoint.z();
         // rotate only in the x-y plane

-        

+

         double x = oldPoint.x();
         double y = oldPoint.y();
         double theta = trackRotationAngle;

-        
-        double xPrime = x*cos(theta) + y*sin(theta);
-        double yPrime = - x*sin(theta) + y*cos(theta);

+
+        double xPrime = x * cos(theta) + y * sin(theta);
+        double yPrime = -x * sin(theta) + y * cos(theta);

         SpacePoint result = new CartesianPoint(xPrime, yPrime, z);
         return result;
     }
 

-    

     public SpacePoint referencePoint() {
         updateRefPoint();
         return _referencePoint;

@@ -384,7 +442,7 @@

         double px = t.getPX();
         double py = t.getPY();
         double pt = sqrt(px * px + py * py);

-        

+

         _parameters = new ZvParameters();
         _parameters.kappa = 1 / pt;
         // FIXME hardcoded order - dangerous

@@ -396,50 +454,211 @@

         _parameters.z = t.getReferencePointZ();
         updateMomentum();
         updateRefPoint();

-        if (! _refPointIsPCA) {

+        if (!_refPointIsPCA) {

             System.err.println("WARNING ! Reference Point should be PCA");
         }
         if (px < 0.0)

-            trackRotationAngle = acos(py/pt);

+            trackRotationAngle = acos(py / pt);

         else

-            trackRotationAngle = - acos(py/pt);

+            trackRotationAngle = -acos(py / pt);

         return;
     }

-    
-    /** 

+
+    /**

      * Calculates the momentum from the parameters
      */
     private void updateMomentum() {

-        double pt = 1/_parameters.kappa;

+        double pt = 1 / _parameters.kappa;

         double px = pt * Math.cos(_parameters.phi);
         double py = pt * Math.sin(_parameters.phi);
         double pz = pt * _parameters.s;
         _momentum = new BasicHep3Vector(px, py, pz);

-        

+

     }

-    

+

     /**
      * Calculates the reference point from the parameters
      */
     private void updateRefPoint() {
         _referencePoint = new CartesianPoint(_parameters.x, _parameters.y, _parameters.z);
     }

-    

+

     /**
      * Fills the errorMatrix of a ZvTrack from a double[5][5] of a Track

+     * 

      * @param trackMatrix
      */
     // TODO implementation
     private void convert2ZvErrorMatrix(double[][] trackMatrix) {

-        

+

     }

-    

+

     /**
      * converts the ZvErrorMatrix back to a Track-compatible Matrix

+     * 

      * @return the error Matrix in a Track-compatible way
      */
     // TODO implementation
     private double[][] convert2TrackErrorMatrix() {
         return _errorMatrix;
     }

+
+    /**
+     * Calculates the lengthParameter s at which the ZvTrack is closest to the
+     * given SpacePoint
+     * 
+     * @param point
+     *            The point in space to which to calculate the distance
+     * @return the length parameter
+     */
+    // Apologies about the form. Automatic code generation by Maple(C)
+    double closestParabolaPoint(SpacePoint point) {
+        double x1 = point.x();
+        double y1 = point.y();
+        double z1 = point.z();
+        double kappa = _parameters.kappa;
+        double tanLambda = _parameters.s;
+        double theta = trackRotationAngle;
+        double x0 = rotateLocation(_referencePoint).x();
+        double z0 = _referencePoint.z();
+        return ((Math.pow(0.54e2
+                * kappa
+                * y1
+                * Math.cos(theta)
+                - 0.54e2
+                * kappa
+                * tanLambda
+                * z0
+                + 0.54e2
+                * kappa
+                * tanLambda
+                * z1
+                - 0.54e2
+                * kappa
+                * x1
+                * Math.sin(theta)
+                + 0.6e1
+                * Math.sqrt(0.6e1 + 0.72e2 * tanLambda * tanLambda * x0 * kappa + 0.81e2 * kappa
+                        * kappa * tanLambda * tanLambda * z0 * z0 + 0.72e2 * tanLambda * tanLambda
+                        * x0 * x0 * kappa * kappa + 0.36e2 * Math.pow(tanLambda, 0.4e1) * x0
+                        * kappa + 0.36e2 * x0 * kappa + 0.48e2 * Math.pow(x0, 0.3e1)
+                        * Math.pow(kappa, 0.3e1) + 0.72e2 * x0 * x0 * kappa * kappa + 0.6e1
+                        * Math.pow(tanLambda, 0.6e1) + 0.18e2 * Math.pow(tanLambda, 0.4e1) + 0.18e2
+                        * tanLambda * tanLambda + 0.72e2 * y1 * y1 * kappa * kappa + 0.81e2 * kappa
+                        * kappa * x1 * x1 + 0.144e3 * y1 * kappa * kappa * Math.sin(theta) * x1
+                        * Math.cos(theta) * tanLambda * tanLambda + 0.162e3 * kappa * kappa * y1
+                        * Math.cos(theta) * tanLambda * z1 + 0.288e3 * y1 * Math.pow(kappa, 0.3e1)
+                        * Math.sin(theta) * x1 * Math.cos(theta) * x0 - 0.162e3 * kappa * kappa
+                        * y1 * Math.cos(theta) * tanLambda * z0 - 0.144e3 * x1 * kappa * kappa
+                        * Math.cos(theta) * tanLambda * tanLambda * x0 - 0.162e3 * kappa * kappa
+                        * tanLambda * z1 * x1 * Math.sin(theta) + 0.162e3 * kappa * kappa
+                        * tanLambda * z0 * x1 * Math.sin(theta) - 0.144e3 * y1
+                        * Math.pow(kappa, 0.3e1) * Math.sin(theta) * x1 * x1
+                        * Math.pow(Math.cos(theta), 0.2e1) - 0.144e3 * y1 * kappa * kappa
+                        * Math.sin(theta) * tanLambda * tanLambda * x0 - 0.18e2 * y1 * kappa
+                        * kappa * Math.sin(theta) * x1 * Math.cos(theta) - 0.144e3
+                        * Math.pow(Math.cos(theta), 0.2e1) * y1 * y1 * Math.pow(kappa, 0.3e1) * x0
+                        + 0.144e3 * Math.pow(Math.cos(theta), 0.3e1) * y1 * y1
+                        * Math.pow(kappa, 0.3e1) * x1 - 0.144e3 * y1 * y1 * Math.pow(kappa, 0.3e1)
+                        * x1 * Math.cos(theta) - 0.72e2 * Math.pow(Math.cos(theta), 0.2e1) * y1
+                        * y1 * kappa * kappa * tanLambda * tanLambda + 0.48e2 * Math.sin(theta)
+                        * Math.pow(y1, 0.3e1) * Math.pow(kappa, 0.3e1)
+                        * Math.pow(Math.cos(theta), 0.2e1) + 0.72e2 * y1 * y1 * kappa * kappa
+                        * tanLambda * tanLambda + 0.81e2 * kappa * kappa * tanLambda * tanLambda
+                        * z1 * z1 - 0.48e2 * Math.pow(x1, 0.3e1) * Math.pow(kappa, 0.3e1)
+                        * Math.pow(Math.cos(theta), 0.3e1) + 0.9e1 * kappa * kappa * y1 * y1
+                        * Math.pow(Math.cos(theta), 0.2e1) - 0.9e1 * x1 * x1 * kappa * kappa
+                        * Math.pow(Math.cos(theta), 0.2e1) - 0.36e2 * y1 * kappa * Math.sin(theta)
+                        - 0.36e2 * x1 * kappa * Math.cos(theta) + 0.144e3 * y1 * y1
+                        * Math.pow(kappa, 0.3e1) * x0 - 0.48e2 * Math.sin(theta)
+                        * Math.pow(y1, 0.3e1) * Math.pow(kappa, 0.3e1) - 0.144e3 * y1 * kappa
+                        * kappa * Math.sin(theta) * x0 - 0.144e3 * x1 * kappa * kappa
+                        * Math.cos(theta) * x0 - 0.36e2 * y1 * kappa * Math.sin(theta)
+                        * Math.pow(tanLambda, 0.4e1) - 0.72e2 * y1 * kappa * Math.sin(theta)
+                        * tanLambda * tanLambda - 0.144e3 * y1 * Math.pow(kappa, 0.3e1)
+                        * Math.sin(theta) * x0 * x0 - 0.162e3 * kappa * kappa * tanLambda
+                        * tanLambda * z0 * z1 - 0.144e3 * x1 * Math.pow(kappa, 0.3e1)
+                        * Math.cos(theta) * x0 * x0 - 0.36e2 * x1 * kappa * Math.cos(theta)
+                        * Math.pow(tanLambda, 0.4e1) + 0.144e3 * x1 * x1 * Math.pow(kappa, 0.3e1)
+                        * Math.pow(Math.cos(theta), 0.2e1) * x0 + 0.72e2 * x1 * x1 * kappa * kappa
+                        * Math.pow(Math.cos(theta), 0.2e1) * tanLambda * tanLambda - 0.72e2 * x1
+                        * kappa * Math.cos(theta) * tanLambda * tanLambda), 0.2e1 / 0.3e1)
+                + 0.12e2
+                * y1
+                * kappa
+                * Math.sin(theta)
+                + 0.12e2
+                * x1
+                * kappa
+                * Math.cos(theta)
+                - 0.6e1 * tanLambda * tanLambda - 0.12e2 * x0 * kappa - 0.6e1)
+                * Math.pow(0.54e2
+                        * kappa
+                        * y1
+                        * Math.cos(theta)
+                        - 0.54e2
+                        * kappa
+                        * tanLambda
+                        * z0
+                        + 0.54e2
+                        * kappa
+                        * tanLambda
+                        * z1
+                        - 0.54e2
+                        * kappa
+                        * x1
+                        * Math.sin(theta)
+                        + 0.6e1
+                        * Math.sqrt(0.6e1 + 0.72e2 * tanLambda * tanLambda * x0 * kappa + 0.81e2
+                                * kappa * kappa * tanLambda * tanLambda * z0 * z0 + 0.72e2
+                                * tanLambda * tanLambda * x0 * x0 * kappa * kappa + 0.36e2
+                                * Math.pow(tanLambda, 0.4e1) * x0 * kappa + 0.36e2 * x0 * kappa
+                                + 0.48e2 * Math.pow(x0, 0.3e1) * Math.pow(kappa, 0.3e1) + 0.72e2
+                                * x0 * x0 * kappa * kappa + 0.6e1 * Math.pow(tanLambda, 0.6e1)
+                                + 0.18e2 * Math.pow(tanLambda, 0.4e1) + 0.18e2 * tanLambda
+                                * tanLambda + 0.72e2 * y1 * y1 * kappa * kappa + 0.81e2 * kappa
+                                * kappa * x1 * x1 + 0.144e3 * y1 * kappa * kappa * Math.sin(theta)
+                                * x1 * Math.cos(theta) * tanLambda * tanLambda + 0.162e3 * kappa
+                                * kappa * y1 * Math.cos(theta) * tanLambda * z1 + 0.288e3 * y1
+                                * Math.pow(kappa, 0.3e1) * Math.sin(theta) * x1 * Math.cos(theta)
+                                * x0 - 0.162e3 * kappa * kappa * y1 * Math.cos(theta) * tanLambda
+                                * z0 - 0.144e3 * x1 * kappa * kappa * Math.cos(theta) * tanLambda
+                                * tanLambda * x0 - 0.162e3 * kappa * kappa * tanLambda * z1 * x1
+                                * Math.sin(theta) + 0.162e3 * kappa * kappa * tanLambda * z0 * x1
+                                * Math.sin(theta) - 0.144e3 * y1 * Math.pow(kappa, 0.3e1)
+                                * Math.sin(theta) * x1 * x1 * Math.pow(Math.cos(theta), 0.2e1)
+                                - 0.144e3 * y1 * kappa * kappa * Math.sin(theta) * tanLambda
+                                * tanLambda * x0 - 0.18e2 * y1 * kappa * kappa * Math.sin(theta)
+                                * x1 * Math.cos(theta) - 0.144e3 * Math.pow(Math.cos(theta), 0.2e1)
+                                * y1 * y1 * Math.pow(kappa, 0.3e1) * x0 + 0.144e3
+                                * Math.pow(Math.cos(theta), 0.3e1) * y1 * y1
+                                * Math.pow(kappa, 0.3e1) * x1 - 0.144e3 * y1 * y1
+                                * Math.pow(kappa, 0.3e1) * x1 * Math.cos(theta) - 0.72e2
+                                * Math.pow(Math.cos(theta), 0.2e1) * y1 * y1 * kappa * kappa
+                                * tanLambda * tanLambda + 0.48e2 * Math.sin(theta)
+                                * Math.pow(y1, 0.3e1) * Math.pow(kappa, 0.3e1)
+                                * Math.pow(Math.cos(theta), 0.2e1) + 0.72e2 * y1 * y1 * kappa
+                                * kappa * tanLambda * tanLambda + 0.81e2 * kappa * kappa
+                                * tanLambda * tanLambda * z1 * z1 - 0.48e2 * Math.pow(x1, 0.3e1)
+                                * Math.pow(kappa, 0.3e1) * Math.pow(Math.cos(theta), 0.3e1) + 0.9e1
+                                * kappa * kappa * y1 * y1 * Math.pow(Math.cos(theta), 0.2e1)
+                                - 0.9e1 * x1 * x1 * kappa * kappa
+                                * Math.pow(Math.cos(theta), 0.2e1) - 0.36e2 * y1 * kappa
+                                * Math.sin(theta) - 0.36e2 * x1 * kappa * Math.cos(theta) + 0.144e3
+                                * y1 * y1 * Math.pow(kappa, 0.3e1) * x0 - 0.48e2 * Math.sin(theta)
+                                * Math.pow(y1, 0.3e1) * Math.pow(kappa, 0.3e1) - 0.144e3 * y1
+                                * kappa * kappa * Math.sin(theta) * x0 - 0.144e3 * x1 * kappa
+                                * kappa * Math.cos(theta) * x0 - 0.36e2 * y1 * kappa
+                                * Math.sin(theta) * Math.pow(tanLambda, 0.4e1) - 0.72e2 * y1
+                                * kappa * Math.sin(theta) * tanLambda * tanLambda - 0.144e3 * y1
+                                * Math.pow(kappa, 0.3e1) * Math.sin(theta) * x0 * x0 - 0.162e3
+                                * kappa * kappa * tanLambda * tanLambda * z0 * z1 - 0.144e3 * x1
+                                * Math.pow(kappa, 0.3e1) * Math.cos(theta) * x0 * x0 - 0.36e2 * x1
+                                * kappa * Math.cos(theta) * Math.pow(tanLambda, 0.4e1) + 0.144e3
+                                * x1 * x1 * Math.pow(kappa, 0.3e1)
+                                * Math.pow(Math.cos(theta), 0.2e1) * x0 + 0.72e2 * x1 * x1 * kappa
+                                * kappa * Math.pow(Math.cos(theta), 0.2e1) * tanLambda * tanLambda
+                                - 0.72e2 * x1 * kappa * Math.cos(theta) * tanLambda * tanLambda),
+                        -0.1e1 / 0.3e1) / kappa / 0.6e1);
+    }

 }
lcsim/src/org/lcsim/recon/vertexing/zvtop4
ZvTop.java 1.9 -> 1.10
diff -u -r1.9 -r1.10
--- ZvTop.java	27 Jul 2005 03:03:08 -0000	1.9
+++ ZvTop.java	1 Aug 2005 18:09:50 -0000	1.10
@@ -36,6 +36,8 @@

  * @author jstrube
  */
 public class ZvTop {

+    // stores the f_i(r)*f_j(r), which are later translated into the
+    // V(r_ij)

     ZvMaximumMatrix maximumMatrix;
     List<ZvTrack> trackList;
     private ZvFitter vtxFitter;

@@ -134,14 +136,15 @@

         return;
     }
       

-    // TODO profiling    

+    /**
+     * Finds the Maxima
+     * @param list
+     * @return
+     */

     public List<ZvVertex> findVertices(List<Track> list) {
         trackList.clear();
         System.err.println("entering findVertives\n");

-        // empty track to represent IP
-        // FIXME Put in IP as constraint
-//        ZvTrack ipTrack = new ZvTrack();
-//        trackList.add(ipTrack);

+//         FIXME Put in IP as constraint. At the moment, no vertex contains IP        

         for (Track iTrack : list) {
             trackList.add(new ZvTrack(iTrack));
         }

@@ -157,13 +160,19 @@

             }
         }
         assignResolution();

-        transformMaxima();

+        findGlobalMaxima();

         clusterCandidates();
         System.err.println("done findVertices\n");
         return makeVertices();
     }
       

-    // units are mm

+    /**
+     * Calculates the global overlap at a given point in space
+     * The global overlap is the sum S of all track tubes minus
+     * the sum of the squares of the track tubes divided by S 
+     * @param location The point in space at which to evaluate the function
+     * @return The value at that point 
+     */

     public double overlap(SpacePoint location) {
         double numerator = 0;
         double denominator = 0;

@@ -243,10 +252,115 @@

         }
     }
     

-    // TODO dummy function
-    // Supposed to 'transform the r_ij to the nearest maximum in V(r)'
-    // Does transformation change the order in the list ?
-    private void transformMaxima() {

+    /**
+     * Transforms each of the maxima in the two-track overlap 
+     * to the nearest global maximum in V(r).
+     * Finds the global maxima first.
+     */
+    // TODO maybe the maxima don't even have to be transformed
+    // TODO the sequence of for-loops is a little clumsy
+    private void findGlobalMaxima() {
+        // search in micrometer steps in xy
+        final double xyStep = 0.001;
+        // search in 10 micrometer steps in z
+        final double zStep = 0.01;
+        final int maxSteps = 500;
+        for (ZvMaximum iMax : maximumMatrix.getMaxima()) {
+            SpacePoint searchPoint = iMax.getLocation();
+            SpacePoint maxPoint = new SpacePoint(searchPoint);
+            double xPos = searchPoint.x();
+            double yPos = searchPoint.y();
+            double zPos = searchPoint.z();
+            double globalMax = iMax.getValue();
+            // perform the iteration __consecutively__, not simultaneously !
+            // for every direction: go in positive direction,
+            // reset to initial position and go in negative direction
+            // +x
+            for (int x=1; x<=maxSteps; ++x) {
+                xPos += xyStep;
+                searchPoint = new CartesianPoint(xPos, yPos, zPos);
+                double overlap = overlap(searchPoint);
+                if (overlap > globalMax) {
+                    globalMax = overlap;
+                    maxPoint = new SpacePoint(searchPoint);
+                } else {
+                    xPos -= xyStep;
+                    break;
+                }
+            }
+            xPos = iMax.getLocation().x();
+            // -x
+            for (int x=1; x<=maxSteps; ++x) {
+                xPos -= xyStep;
+                searchPoint = new CartesianPoint(xPos, yPos, zPos);
+                double overlap = overlap(searchPoint);
+                if (overlap > globalMax) {
+                    globalMax = overlap;
+                    maxPoint = new SpacePoint(searchPoint);
+                } else {
+                    xPos += xyStep;
+                    break;
+                }
+            }
+            xPos = maxPoint.x();
+            // +y
+            for (int y=1; y<=maxSteps; ++y) {
+                yPos += xyStep;
+                searchPoint = new CartesianPoint(xPos, yPos, zPos);
+                double overlap = overlap(searchPoint);
+                if (overlap > globalMax) {
+                    globalMax = overlap;
+                    maxPoint = new SpacePoint(searchPoint);
+                } else {
+                    yPos -= xyStep;
+                    break;
+                }
+            }
+            yPos = iMax.getLocation().y();
+            // -y
+            for (int y=1; y<=maxSteps; ++y) {
+                yPos -= xyStep;
+                searchPoint = new CartesianPoint(xPos, yPos, zPos);
+                double overlap = overlap(searchPoint);
+                if (overlap > globalMax) {
+                    globalMax = overlap;
+                    maxPoint = new SpacePoint(searchPoint);
+                } else {
+                    yPos += xyStep;
+                    break;
+                }
+            }
+            yPos = maxPoint.y();
+            // +z
+            for (int z=1; z<=maxSteps; ++z) {
+                zPos += zStep;
+                searchPoint = new CartesianPoint(xPos, yPos, zPos);
+                double overlap = overlap(searchPoint);
+                if (overlap > globalMax) {
+                    globalMax = overlap;
+                    maxPoint = new SpacePoint(searchPoint);
+                } else {
+                    zPos -= zStep;
+                    break;
+                }
+            }
+            zPos = iMax.getLocation().z();
+            // -z
+            for (int z=1; z<=maxSteps; ++z) {
+                zPos -= zStep;
+                searchPoint = new CartesianPoint(xPos, yPos, zPos);
+                double overlap = overlap(searchPoint);
+                if (overlap > globalMax) {
+                    globalMax = overlap;
+                    maxPoint = new SpacePoint(searchPoint);
+                } else {
+                    zPos += zStep;
+                    break;
+                }
+            }
+            iMax.setLocation(maxPoint);
+        }
+        

     }
     
     /**

@@ -282,7 +396,6 @@

      * always assigned to the vertex with the highest value of V(r), if possible
      * @return 
      */

-    // TODO

     private List<ZvVertex> makeVertices() {
         List<ZvVertex> l = new ArrayList<ZvVertex>();
         Set<ZvTrack> unavailableTracks = new HashSet<ZvTrack>();
lcsim/src/org/lcsim/recon/vertexing/zvtop4
ZvMaximum.java 1.1 -> 1.2
diff -u -r1.1 -r1.2
--- ZvMaximum.java	19 Apr 2005 22:22:56 -0000	1.1
+++ ZvMaximum.java	1 Aug 2005 18:09:50 -0000	1.2
@@ -25,9 +25,14 @@

         _value = overlap;
     }
     

-    // Set new location
-    // Assume that overlap doesn't change enough for the order of ZvMaxima to change
-    public void setLocation(SpacePoint loc) {

+    /**
+     * Sets the new location of the ZvMaximum (used by clusterCandidates)
+     * @param loc new location of the maximum. It is assumed that the 
+     * overlap doesn't change enough between the two points for the order to
+     * change in the ZvMaximumMatrix
+     * Not intended for public consumption
+     */
+    void setLocation(SpacePoint loc) {

         _location = loc;
     }
CVSspam 0.2.8