lcsim/src/org/lcsim/contrib/JanStrube/zvtop
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
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
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
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
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;
+ }
+}