Print

Print


Commit in lcsim/src/org/lcsim/contrib/JanStrube/zvtop on MAIN
ZvTrack.java+24added 1.1
ZvMaximum.java+67added 1.1
ZvMaximumMatrix.java+72added 1.1
ZvUtil.java+129added 1.1
ZvTop.java+489added 1.1
+781
5 added files
testbed for integration of Kalman fitter to zvtop

lcsim/src/org/lcsim/contrib/JanStrube/zvtop
ZvTrack.java added at 1.1
diff -N ZvTrack.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ZvTrack.java	9 Oct 2006 17:25:27 -0000	1.1
@@ -0,0 +1,24 @@
+/**
+ * @version $Id: ZvTrack.java,v 1.1 2006/10/09 17:25:27 jstrube Exp $
+ */
+package org.lcsim.contrib.JanStrube.zvtop;
+
+import org.lcsim.contrib.JanStrube.tracking.NewTrack;
+import org.lcsim.contrib.JanStrube.tracking.Track;
+import org.lcsim.spacegeom.SpacePoint;
+
+/**
+ * @author jstrube
+ *
+ */
+public class ZvTrack {
+    /**
+     * returns the probability that a track goes through a certain point
+     * Assuming gaussian errors.
+     * After Lyons, p. 47
+     */
+    public static double getProbability(Track t, SpacePoint l) {
+        throw new RuntimeException("Not implemented");
+//        return 0;
+    }
+}

lcsim/src/org/lcsim/contrib/JanStrube/zvtop
ZvMaximum.java added at 1.1
diff -N ZvMaximum.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ZvMaximum.java	9 Oct 2006 17:25:27 -0000	1.1
@@ -0,0 +1,67 @@
+package org.lcsim.contrib.JanStrube.zvtop;
+
+import org.lcsim.contrib.JanStrube.tracking.Track;
+import org.lcsim.spacegeom.SpacePoint;
+
+/**
+ * Maximum in the overlap of two tracks
+ * @author jstrube
+ * @version $Id: ZvMaximum.java,v 1.1 2006/10/09 17:25:27 jstrube Exp $
+ */
+class ZvMaximum implements Comparable {
+    private Track _track_i;
+    private Track _track_j;
+    private SpacePoint _location;
+    private double _value;
+    
+    public ZvMaximum(Track track_i, Track track_j, SpacePoint loc, double overlap) {
+        _track_i = track_i;
+        _track_j = track_j;
+        _location = loc;
+        // TODO for debugging only
+        _value = overlap;
+    }
+    
+    /**
+     * 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;
+    }
+    
+    public SpacePoint getLocation() {
+        return _location;
+    }
+    
+    public double getValue() {
+        return _value;
+    }
+    
+    public Track getTrack_i() {
+        return _track_i;
+    }
+    
+    public Track getTrack_j() {
+        return _track_j;
+    }
+    
+    public int compareTo(Object o) {
+        if (o.getClass() != ZvMaximum.class)
+            throw new IllegalArgumentException("Can only compare ZvMaximum Objects");
+        ZvMaximum max2 = (ZvMaximum) o;
+        if (_value > max2._value)
+            return 1;
+        if (_value < max2._value)
+            return -1;
+        return 0;
+    }
+    
+    public String toString() {
+        String s = String.format("value: %f\nlocation: %s\n", _value, _location);
+        return s;
+    }
+}

lcsim/src/org/lcsim/contrib/JanStrube/zvtop
ZvMaximumMatrix.java added at 1.1
diff -N ZvMaximumMatrix.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ZvMaximumMatrix.java	9 Oct 2006 17:25:27 -0000	1.1
@@ -0,0 +1,72 @@
+package org.lcsim.contrib.JanStrube.zvtop;
+
+import java.util.Comparator;
+import java.util.Iterator;
+import java.util.Set;
+import java.util.SortedMap;
+import java.util.TreeMap;
+import java.util.Map.Entry;
+
+/**
+ * Container of ZvMaxima
+ * This enables for the ZvMaxima to be sorted by V(r)
+ * while they still can be accessed by index (i, j) 
+ * @author jstrube
+ * @version $Id: ZvMaximumMatrix.java,v 1.1 2006/10/09 17:25:27 jstrube Exp $
+ *
+ */
+class ZvMaximumMatrix implements Iterable {
+
+    private class MatrixIndex {
+        MatrixIndex(int i, int j) {
+            _i = i;
+            _j = j;
+        }
+        int _i;
+        int _j;
+    }
+    SortedMap<ZvMaximum, MatrixIndex> _table; 
+       
+    /**
+     * 
+     */
+    ZvMaximumMatrix(Comparator comp) {
+        _table = new TreeMap<ZvMaximum, MatrixIndex>(comp);
+    }
+    
+    void add(ZvMaximum max, int i, int j) {
+        _table.put(max, new MatrixIndex(i, j));
+    }
+
+    /* (non-Javadoc)
+     * @see java.lang.Iterable#iterator()
+     */
+    public Iterator iterator() {
+        return _table.keySet().iterator();
+    }
+    
+    /**
+     * Accessor of ZvMaximum by indices of the constituent tracks in the list
+     * @param i index of the first Track
+     * @param j index of the second Track
+     * @return the ZvMaximum
+     */
+    ZvMaximum get(int i, int j) throws IllegalAccessException {
+        // Slow Access - OK, because primary access is by order in V(r)
+        Set<Entry<ZvMaximum, MatrixIndex>> entries = _table.entrySet();
+        for (Entry<ZvMaximum, MatrixIndex> entry  : entries) {
+            MatrixIndex index = entry.getValue();
+            // Symmetry in indices
+            if ((index._i == i && index._j == j)
+                    || (index._i == j && index._j == i) ) {
+                return entry.getKey();
+            }
+        }
+        throw new IllegalAccessException("ZvMaximumMatrix.get: index not found");
+    }
+    
+    Set<ZvMaximum> getMaxima() {
+        return _table.keySet();
+    }
+
+}

lcsim/src/org/lcsim/contrib/JanStrube/zvtop
ZvUtil.java added at 1.1
diff -N ZvUtil.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ZvUtil.java	9 Oct 2006 17:25:27 -0000	1.1
@@ -0,0 +1,129 @@
+package org.lcsim.contrib.JanStrube.zvtop;
+
+import java.util.Collection;
+
+import org.lcsim.contrib.JanStrube.tracking.HelixSwimmer;
+import org.lcsim.contrib.JanStrube.tracking.Track;
+import org.lcsim.geometry.Detector;
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpacePointVector;
+import org.lcsim.spacegeom.SpaceVector;
+
+import Jama.Matrix;
+
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.subtract;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.dot;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.add;
+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
+ * @version $Id: ZvUtil.java,v 1.1 2006/10/09 17:25:27 jstrube Exp $
+ */
+final class ZvUtil {
+
+    /**
+     * Private constructor ensures the class is never instantiated
+     */
+    private ZvUtil() {
+    }
+
+    /**
+     * Calculates the location where the product of the two "track probabilities"
+     * is at a maximum
+     * @param iTrack
+     * @param jTrack
+     * @return the location of the maximum
+     */
+    public static SpacePoint twoTrackMax(Track iTrack, Track jTrack) {
+        return swimmerTwoTrackMax(iTrack, jTrack);
+    }
+    
+    /*
+     * calculates the maximum of the two track overlap using the swimmer
+     * The error is computed using FIXME
+     */
+    private static SpacePoint swimmerTwoTrackMax(Track iTrack, Track jTrack) {
+        HelixSwimmer swimmer = new HelixSwimmer(5);
+        swimmer.setTrack(iTrack);
+        swimmer.getTrackLengthToPoint(jTrack.getPosition());
+        throw new RuntimeException("Implement me !");
+//        return new SpacePoint();
+    }
+    
+    
+    /**
+     * Sums the momenta of all tracks in the list
+     * @param list
+     * @return the vector sum of the momenta
+     */
+    public static SpaceVector sumTrackMomenta(Collection<Track> list) {
+        double x = 0.;
+        double y = 0.;
+        double z = 0.;
+        for (Track iTrack: list) {
+            double[] p3 = iTrack.getMomentum().v();
+            x += p3[0];
+            y += p3[1];
+            z += p3[2];
+        }
+        return new CartesianVector(x, y, z);
+    }
+    
+    /**
+     * Sums the charge of all tracks in the list
+     * @param list
+     * @return the sum of track charges
+     */
+    public static double sumTrackCharge(Collection<Track> list) {
+        double result = 0;
+        for (Track iTrack : list) {
+            result += iTrack.getCharge();
+        }
+        return result;
+    }
+
+    
+    /**
+     * Projection of error matrix along specific direction <br>
+     * (Does normalization of direction vector.)
+     * @param  A          error matrix to be projected (3x3)
+     * @param  direction  direction to project along (3-dim)
+     * @return square of sigma along direction
+     */
+    /*
+     * TODO Legacy code
+     */
+    public static double projectSigma2(double[][]A, SpaceVector direction)
+    {
+        SpaceVector dir = unit(direction);
+        double[] d = new double[] {dir.x(), dir.y(), dir.z()};
+        double sig2 = 0;
+        for (int i=0; i<3; i++) {
+            for (int j=0; j<3; j++) 
+                sig2 += d[i]*A[i][j]*d[j];
+        }
+        return sig2;
+    }
+           
+    /**
+     * calculate p/tangent unit vector
+     */
+    public static SpaceVector getUnitTangent(double[] hlxPar)
+    {
+        int tanLambda = Track.ParameterName.tanLambda.ordinal();
+        int phi = Track.ParameterName.phi0.ordinal();
+        double norm = Math.sqrt(1.+hlxPar[tanLambda]*hlxPar[tanLambda]);
+        double x = Math.cos(hlxPar[phi])/norm;
+        double y = Math.sin(hlxPar[phi])/norm;
+        double z = hlxPar[tanLambda]/norm;
+        return new CartesianVector(x, y, z);
+    }
+    
+}

lcsim/src/org/lcsim/contrib/JanStrube/zvtop
ZvTop.java added at 1.1
diff -N ZvTop.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ZvTop.java	9 Oct 2006 17:25:27 -0000	1.1
@@ -0,0 +1,489 @@
+package org.lcsim.contrib.JanStrube.zvtop;
+
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.add;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.dot;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.multiply;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.subtract;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.unit;
+import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.angle;
+import static org.lcsim.contrib.JanStrube.zvtop.ZvUtil.twoTrackMax;
+
+import java.util.ArrayList;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import java.util.SortedSet;
+import java.util.TreeSet;
+
+import org.lcsim.contrib.JanStrube.tracking.Track;
+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 org.lcsim.spacegeom.SpacePointVector;
+import org.lcsim.spacegeom.SpaceVector;
+
+// TODO think about implementation in terms of static members
+/**
+ * Main class to handle the vertexing Vertices are found by the function
+ * @see #findVertices(List)
+ * @author jstrube
+ * @version $Id: ZvTop.java,v 1.1 2006/10/09 17:25:27 jstrube Exp $
+ */
+public class ZvTop {
+    /**
+     * @author jstrube
+     */
+    class MaxComp implements Comparator {
+
+        /**
+         * 
+         */
+        MaxComp() {
+            super();
+            // TODO Auto-generated constructor stub
+        }
+
+        /*
+         * (non-Javadoc)
+         * @see java.util.Comparator#compare(T, T)
+         */
+        public int compare(Object arg1, Object arg2) {
+            if (arg1.getClass() != ZvMaximum.class
+                    && arg2.getClass() != ZvMaximum.class)
+                throw new IllegalArgumentException(
+                        "Can only compare ZvMaximum Objects");
+            ZvMaximum max1 = (ZvMaximum) arg1;
+            ZvMaximum max2 = (ZvMaximum) arg2;
+            double overlap1 = max1.getValue();
+            double overlap2 = max2.getValue();
+            if (overlap1 > overlap2)
+                return 1;
+            if (overlap1 < overlap2)
+                return -1;
+            return 0;
+        }
+
+    }
+
+    // if a two-track overlap is less than this cut, do not add max to the list
+    private double _globalOverlapCut = 1E-10;
+    private SpaceVector _jetAxis;
+    private double _maxDistFromIP = 25;
+    // List of vertex candidates.
+    // The two-track maxima are transformed to the nearest location of a maximum
+    // in the overlap function.
+    // The highest of those maxima is taken as the location of the vertex and
+    // all tracks belonging to ZvMaxima that are unresolved from this highest
+    // maximum are added to the vertex.
+    // At this stage tracks can be assigned to more than one vertex.
+    // The fitVertices function takes care of assigning tracks uniquely to a
+    // vertex.
+    // TODO need to distinguish between pruning and fitting
+    private List<Vertex> _vertexCandidateList;
+    // Special weight for the angle of the spatial point (argument of overlap)
+    // wrt. the jet axis
+    private double angularWeight = 5;
+    private double chiSquareCut = 15;
+    // These data members must be final, because changing their values
+    // has an effect on many of the other classes
+    // Special weight for the IP in the overlap function
+    private double ipWeight = 1;
+    // resolution function is expensive in CPU
+    private Map<ZvMaximum, Set<ZvMaximum>> isResolvedFromMap;
+
+    // Number of points to check along the line from point 1 to point 2
+    int iterationMax = 3;
+
+    // stores the f_i(r)*f_j(r), which are later translated into the
+    // V(r_ij)
+    ZvMaximumMatrix maximumMatrix;
+    private double resolvedCut = 0.7;
+
+    private double sigmaL;
+
+    private double sigmaT;
+
+    List<Track> trackList;
+
+    private Map<ZvMaximum, Set<ZvMaximum>> unresolvedMaximaMap;
+
+    private Fitter vtxFitter;
+
+    /**
+     * Constructor Creates a list of ZvTracks from the argument list. Calculates
+     * the overlap and resolution functions for the tracks
+     * @param jetAxisValue
+     */
+    public ZvTop(SpaceVector jetAxisValue, Fitter f) {
+        // FIXME Using arbitrary values for the measurement errors
+        sigmaT = 0.01;
+        sigmaL = 0.1;
+        _jetAxis = jetAxisValue;
+        maximumMatrix = new ZvMaximumMatrix(new MaxComp());
+        isResolvedFromMap = new HashMap<ZvMaximum, Set<ZvMaximum>>();
+        unresolvedMaximaMap = new HashMap<ZvMaximum, Set<ZvMaximum>>();
+        trackList = new ArrayList<Track>();
+        _vertexCandidateList = new ArrayList<Vertex>();
+        vtxFitter = f;
+        return;
+    }
+
+    /**
+     * Calculates the resolution for two points in space. Two points in space
+     * are resolved, if the maximum overlap along the line that connects them
+     * divided by the greater of the overlap values at the two points is less
+     * than a cutoff value R0
+     * @param sp1 point 1
+     * @param sp2 point 2
+     * @return true, if the resolution function is less than the cutoff, false
+     *         otherwise
+     */
+    public boolean areResolved(SpacePoint sp1, SpacePoint sp2) {
+        double overlap1 = overlap(sp1);
+        double overlap2 = overlap(sp2);
+        double minOverlap = 200;
+        SpaceVector vec = subtract(sp2, sp1);
+        for (int i = 1; i < iterationMax; ++i) {
+            SpacePoint rVec = add(sp1, multiply(vec, 1.0*i / iterationMax));
+            double overlap = overlap(rVec);
+            if (overlap < minOverlap)
+                minOverlap = overlap;
+        }
+        if (overlap1 < minOverlap)
+            minOverlap = overlap1;
+        if (overlap2 < minOverlap)
+            minOverlap = overlap2;
+//        System.out.println("Value: " + minOverlap
+//                / ((overlap2 < overlap1) ? overlap2 : overlap1));
+        return minOverlap / ((overlap2 < overlap1) ? overlap2 : overlap1) < resolvedCut;
+    }
+
+    /**
+     * Calculates for each ZvMaximum whether of not it is resolved from other
+     * ZvMaxima. Fills two maps to cache those values. Each ZvMaximum is mapped
+     * to a set of resolved/unresolved ZvMaxima
+     */
+    private void assignResolution() {
+//        System.out.println("assignResolution:");
+        final Set<ZvMaximum> maxSet = maximumMatrix.getMaxima();
+        for (ZvMaximum iMax : maxSet) {
+            Set<ZvMaximum> isResolvedFromImax = new HashSet<ZvMaximum>();
+            Set<ZvMaximum> isNotResolvedFromImax = new HashSet<ZvMaximum>();
+            for (ZvMaximum jMax : maxSet) {
+                if (iMax == jMax)
+                    continue;
+//                System.out.printf("%s and %s ", iMax, jMax);
+                if (areResolved(iMax.getLocation(), jMax.getLocation())) {
+                    isResolvedFromImax.add(jMax);
+//                    System.out.println("are resolved");
+                } else {
+                    isNotResolvedFromImax.add(jMax);
+//                    System.out.println("are not resolved");
+                }
+            }
+            isResolvedFromMap.put(iMax, isResolvedFromImax);
+            unresolvedMaximaMap.put(iMax, isNotResolvedFromImax);
+        }
+//        System.out.printf("isResolvedFromMap: %d\tunresolvedMaximaMap: %d\n",
+//                isResolvedFromMap.size(), unresolvedMaximaMap.size());
+//        System.out.println("ZvMaxima");
+//        for (ZvMaximum iMax : unresolvedMaximaMap.keySet()) {
+//            System.out.println(iMax);
+//            System.out.printf("Is resolved from %d other Maxima\n",
+//                    isResolvedFromMap.get(iMax).size());
+//            System.out.printf("Is unresolved from %d other Maxima\n",
+//                    unresolvedMaximaMap.get(iMax).size());
+//        }
+//        System.out.println("done assignResolution");
+    }
+
+    /**
+     * Creates a list of Vertex candidates. Each candidate consists of a number
+     * of track-IP overlaps that are not spatially resolved. The assignment is
+     * unique in the sense that each overlap is assigned to exactly one
+     * candidate.
+     */
+    private void clusterCandidates() {
+//        System.err.println("entering clusterCandidates");
+        SortedSet<ZvMaximum> isAvailable = new TreeSet<ZvMaximum>(maximumMatrix
+                .getMaxima());
+        // each cluster consists of several ZvMaxima
+        // add the maximum and the set of unresolved maxima to the list of
+        // candidates
+        // __is__ already sorted by overlap value
+        while (!isAvailable.isEmpty()) {
+            ZvMaximum highestRemaining = isAvailable.first();
+            Set<ZvMaximum> unresolvedSet = unresolvedMaximaMap
+                    .get(highestRemaining);
+            unresolvedSet.add(highestRemaining);
+            Vertex vtx = new Vertex();
+            Set<Track> tracks = new HashSet<Track>();
+            for (ZvMaximum imax : unresolvedSet) {
+                tracks.add(imax.getTrack_i());
+                tracks.add(imax.getTrack_j());
+            }
+            _vertexCandidateList.add(vtx);
+            isAvailable.removeAll(unresolvedSet);
+        }
+//        System.err.printf("done clusterCandidates: %d Candidates",
+//                _vertexCandidateList.size());
+        return;
+    }
+
+    /**
+     * 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() {
+//        System.out.println("Entering 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 until maximum is
+            // found,
+            // reset to initial position and go in negative direction
+            // +x
+            for (int x = 1; x <= maxSteps; ++x) {
+                xPos += xyStep;
+                // TODO this is maybe where a double[] is acceptable
+                searchPoint = new CartesianPoint(xPos, yPos, zPos);
+                double overlap = overlap(searchPoint);
+                if (overlap > globalMax) {
+                    globalMax = overlap;
+                    maxPoint = new SpacePoint(searchPoint);
+                } else {
+                    xPos -= xyStep;
+                    break;
+                }
+            }
+            // reset to original pos
+            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;
+                }
+            }
+            // search in x is done. set to found maximum
+            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;
+                }
+            }
+            // reset to original pos
+            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;
+                }
+            }
+            // search in y is done. set to found maximum
+            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;
+                }
+            }
+            // reset to original pos
+            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;
+                }
+            }
+            // search in x,y,z is done. set to found maximum
+            iMax.setLocation(maxPoint);
+        }
+//        System.out.println("Global Maxima done");
+    }
+
+    /**
+     * Finds the Maxima
+     * @param list List of tracks
+     * @return The list of vertices
+     */
+    public List<Vertex> findVertices(List<Track> list) {
+        trackList.clear();
+//        System.err
+//                .printf("entering findVertives with %d Tracks\n", list.size());
+        // FIXME Put in IP as constraint. At the moment, no vertex contains IP
+        for (Track iTrack : list) {
+            trackList.add(iTrack);
+        }
+        int size = trackList.size();
+        for (int i = 0; i < size; ++i) {
+            Track iTrack = trackList.get(i);
+            for (int j = i + 1; j < size; ++j) {
+                Track jTrack = trackList.get(j);
+                SpacePoint location = twoTrackMax(iTrack, jTrack);
+                double overlap = overlap(location);
+                // consider only significant overlap
+                if (overlap < _globalOverlapCut)
+                    continue;
+                ZvMaximum max = new ZvMaximum(iTrack, jTrack, location, overlap);
+                maximumMatrix.add(max, i, j);
+            }
+        }
+//        System.out.println("maximumMatrix: " + maximumMatrix._table.size());
+        findGlobalMaxima();
+        assignResolution();
+        clusterCandidates();
+//        System.err.println("done findVertices\n");
+        return fitVertices();
+    }
+
+    /**
+     * fits the tracks to the vertices. A ZvVertex can only have tracks with a
+     * chi2 contribution below a certain cut value. They are assigned from a
+     * sorted set of maxima, so that the tracks are always assigned to the
+     * vertex with the highest value of V(r), if possible
+     * @return A List of Vertices
+     */
+    private List<Vertex> fitVertices() {
+        List<Vertex> l = new ArrayList<Vertex>();
+        Set<Track> unavailableTracks = new HashSet<Track>();
+//        System.err.printf("found %d Candidates\n", _vertexCandidateList.size());
+        for (Vertex iCand : _vertexCandidateList) {
+//            System.err.printf(iCand.toString());
+        }
+        // The Maxima are clustered according to the resolution function
+        // A Cluster is a Set of ZvMaxima
+        // Add all tracks from a cluster to the vertex
+        // The vertex decides based on the chiSquared value
+        // which of those tracks it wants to keep
+        for (Vertex iVtx : _vertexCandidateList) {
+//            System.out.println("using chi2 cut of " + chiSquareCut);
+            iVtx.setChi2(chiSquareCut);
+            for (Track t : unavailableTracks) {
+                iVtx.removeTrack(t);
+            }
+            // FIXME using a vertex as the container isn't so good
+            // because a) tracks can drop out
+            // so that a newly created vertex should be returned...
+            // that makes for messy bookkeeping
+            vtxFitter.fit(iVtx.getTracks(), iVtx.origin());
+            l.add(iVtx);
+            unavailableTracks.addAll(iVtx.getTracks());
+        }
+        return l;
+    }
+
+    public double getGlobalOverlapCut() {
+        return _globalOverlapCut;
+    }
+
+    /**
+     * 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;
+        // calculate the angular weight
+        // the angle is taken not from (000),
+        // but from a point -0.1 mm along the jet axis
+        double[] origin = multiply(unit(_jetAxis), -0.1).v();
+        SpacePointVector vec2Location = new SpacePointVector(
+                new CartesianPoint(origin[0], origin[1], origin[2]), location);
+        double angle = angle(_jetAxis, vec2Location);
+        double angularFactor = Math.exp(-angularWeight * angle * angle);
+        if (dot(unit(_jetAxis), vec2Location.getDirection()) > _maxDistFromIP) {
+            return 0;
+        }
+        // FIXME dealing with IP
+        double ipFactor = 0; // ipWeight *
+        // trackList.get(0).getTubeValue(location);
+        denominator += ipFactor;
+        numerator += 0; // ipFactor * trackList.get(0).getTubeValue(location);
+        for (int i = 0; i < trackList.size(); ++i) {
+            Track iTrack = trackList.get(i);
+            double iTubeVal = ZvTrack.getProbability(iTrack, location);
+            numerator += iTubeVal * iTubeVal;
+            denominator += iTubeVal;
+        }
+        return angularFactor * (denominator - numerator / denominator);
+    }
+
+    public void setAngularWeight(double w) {
+        angularWeight = w;
+    }
+
+    public void setChiSquareCut(double c) {
+        chiSquareCut = c;
+    }
+
+    public void setGlobalOverlapCut(double globalOverlapCut) {
+        _globalOverlapCut = globalOverlapCut;
+    }
+
+    public void setIpWeight(double w) {
+        ipWeight = w;
+    }
+
+    public void setMaxDistFromIP(double d) {
+        _maxDistFromIP = d;
+    }
+
+    public void setResolvedCut(double c) {
+        resolvedCut = c;
+    }
+}
CVSspam 0.2.8