lcsim/src/org/lcsim/recon/vertexing/zvtop4
diff -u -r1.20 -r1.21
--- ZvTop.java 13 Sep 2005 09:32:51 -0000 1.20
+++ ZvTop.java 15 Sep 2005 00:34:00 -0000 1.21
@@ -1,5 +1,5 @@
/**
- *
+ * @version
*/
package org.lcsim.recon.vertexing.zvtop4;
@@ -33,72 +33,12 @@
// 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
*/
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;
- // 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;
- // Special weight for the angle of the spatial point (argument of overlap)
- // wrt. the jet axis
- private double angularWeight = 10;
- private double resolvedCut = 0;
- private double chiSquareCut = 0;
- private double _maxDistFromIP = 25;
- // Number of points to check along the line from point 1 to point 2
- int iterationMax = 3000;
-
- private double sigmaT;
- private double sigmaL;
-
- public void setMaxDistFromIP(double d) {
- _maxDistFromIP = d;
- }
-
- public void setIpWeight(double w) {
- ipWeight = w;
- }
-
- public void setAngularWeight(double w) {
- angularWeight = w;
- }
-
- public void setResolvedCut(double c) {
- resolvedCut = c;
- }
-
- public void setChiSquareCut(double c) {
- chiSquareCut = c;
- }
-
- private Hep3Vector _jetAxis;
- // resolution function is expensive in CPU
- private Map<ZvMaximum, Set<ZvMaximum>> isResolvedFromMap;
- private Map<ZvMaximum, Set<ZvMaximum>> unresolvedMaximaMap;
-
- // 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<ZvVertex> _vertexCandidateList;
-
/**
* @author jstrube
- *
*/
class MaxComp implements Comparator {
@@ -112,12 +52,13 @@
/*
* (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");
+ 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();
@@ -131,10 +72,53 @@
}
+ // if a two-track overlap is less than this cut, do not add max to the list
+ private double _globalOverlapCut = 1E-10;
+ private Hep3Vector _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<ZvVertex> _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 = 300;
+
+ // 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<ZvTrack> trackList;
+
+ private Map<ZvMaximum, Set<ZvMaximum>> unresolvedMaximaMap;
+
+ private ZvFitter vtxFitter;
+
/**
* Constructor Creates a list of ZvTracks from the argument list. Calculates
* the overlap and resolution functions for the tracks
- *
* @param jetAxisValue
*/
public ZvTop(Hep3Vector jetAxisValue, ZvFitter f) {
@@ -152,106 +136,32 @@
}
/**
- * Finds the Maxima
- *
- * @param list
- * @return
- */
- public List<ZvVertex> 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(new ZvTrack(iTrack));
- }
- int size = trackList.size();
- for (int i = 0; i < size; ++i) {
- ZvTrack iTrack = trackList.get(i);
- for (int j = i + 1; j < size; ++j) {
- ZvTrack jTrack = trackList.get(j);
- SpacePoint location = twoTrackMax(iTrack, jTrack);
- double overlap = overlap(location);
- ZvMaximum max = new ZvMaximum(iTrack, jTrack, location, overlap);
- maximumMatrix.add(max, i, j);
- }
- }
- System.out.println("maximumMatrix: " + maximumMatrix._table.size());
- for (ZvMaximum iMax : maximumMatrix._table.keySet()) {
- System.out.println(iMax);
- }
- assignResolution();
- findGlobalMaxima();
- clusterCandidates();
- System.err.println("done findVertices\n");
- return fitVertices();
- }
-
- /**
- * 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) {
- ZvTrack iTrack = trackList.get(i);
- double iTubeVal = iTrack.getTubeValue(location);
- numerator += iTubeVal * iTubeVal;
- denominator += iTubeVal;
- }
- return angularFactor * (denominator - numerator / denominator);
- }
-
- /**
* 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 minOverlap = 2000;
- // TODO conversion SpacePoint Hep3Vector :-(
+ double overlap1 = overlap(sp1);
+ double overlap2 = overlap(sp2);
+ double minOverlap = 200;
Hep3Vector vec = subtract(sp2, sp1);
- Hep3Vector toSp1 = new BasicHep3Vector(sp1.x(), sp1.y(), sp1.z());
for (int i = 1; i < iterationMax; ++i) {
- // TODO fix the casts
- Hep3Vector rVec = add(toSp1, multiply(vec, i / iterationMax));
- SpacePoint r = new CartesianPoint(rVec.x(), rVec.y(), rVec.z());
- double overlap = overlap(r);
+ SpacePoint rVec = add(sp1, multiply(vec, 1.0*i / iterationMax));
+ double overlap = overlap(rVec);
if (overlap < minOverlap)
minOverlap = overlap;
}
- double overlap1 = overlap(sp1);
- double overlap2 = overlap(sp2);
if (overlap1 < minOverlap)
minOverlap = overlap1;
if (overlap2 < minOverlap)
minOverlap = overlap2;
- // overall minimum / min(overlap1, overlap2) < resolvedCut ?
+ System.out.println("Value: " + minOverlap
+ / ((overlap2 < overlap1) ? overlap2 : overlap1));
return minOverlap / ((overlap2 < overlap1) ? overlap2 : overlap1) < resolvedCut;
}
@@ -261,6 +171,7 @@
* 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>();
@@ -268,26 +179,58 @@
for (ZvMaximum jMax : maxSet) {
if (iMax == jMax)
continue;
- if (areResolved(iMax.getLocation(), jMax.getLocation()))
+ System.out.printf("%s and %s ", iMax, jMax);
+ if (areResolved(iMax.getLocation(), jMax.getLocation())) {
isResolvedFromImax.add(jMax);
- else
+ 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.println("assignResolution:");
- System.out.printf("isResolvedFromMap: %d\tunresolvedMaximaMap: %d\n", isResolvedFromMap.size(), unresolvedMaximaMap.size());
+ 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("List of unresolved Maxima (%d):", unresolvedMaximaMap.get(iMax).size());
- for (ZvMaximum iUnresolvedMax : unresolvedMaximaMap.get(iMax)) {
- System.out.println(iUnresolvedMax);
- }
+ 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);
+ ZvVertex vtx = new ZvVertex(highestRemaining.getLocation(),
+ unresolvedSet);
+ _vertexCandidateList.add(vtx);
+ isAvailable.removeAll(unresolvedSet);
}
- System.out.println("");
+ System.err.printf("done clusterCandidates: %d Candidates",
+ _vertexCandidateList.size());
+ return;
}
/**
@@ -297,6 +240,7 @@
// 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
@@ -310,11 +254,13 @@
double zPos = searchPoint.z();
double globalMax = iMax.getValue();
// perform the iteration __consecutively__, not simultaneously !
- // for every direction: go in positive direction,
+ // 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) {
@@ -325,6 +271,7 @@
break;
}
}
+ // reset to original pos
xPos = iMax.getLocation().x();
// -x
for (int x = 1; x <= maxSteps; ++x) {
@@ -339,6 +286,7 @@
break;
}
}
+ // search in x is done. set to found maximum
xPos = maxPoint.x();
// +y
for (int y = 1; y <= maxSteps; ++y) {
@@ -353,6 +301,7 @@
break;
}
}
+ // reset to original pos
yPos = iMax.getLocation().y();
// -y
for (int y = 1; y <= maxSteps; ++y) {
@@ -367,6 +316,7 @@
break;
}
}
+ // search in y is done. set to found maximum
yPos = maxPoint.y();
// +z
for (int z = 1; z <= maxSteps; ++z) {
@@ -381,6 +331,7 @@
break;
}
}
+ // reset to original pos
zPos = iMax.getLocation().z();
// -z
for (int z = 1; z <= maxSteps; ++z) {
@@ -395,35 +346,45 @@
break;
}
}
+ // search in x,y,z is done. set to found maximum
iMax.setLocation(maxPoint);
}
-
+ System.out.println("Global Maxima done");
}
/**
- * 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.
+ * Finds the Maxima
+ * @param list
+ * @return
*/
- private void clusterCandidates() {
- System.err.println("entering clusterCandidates");
- System.out.println(maximumMatrix.getMaxima());
- 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);
- ZvVertex vtx = new ZvVertex(highestRemaining.getLocation(), unresolvedSet);
- _vertexCandidateList.add(vtx);
- isAvailable.removeAll(unresolvedSet);
+ public List<ZvVertex> 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(new ZvTrack(iTrack));
}
- System.err.printf("done clusterCandidates: %d Candidates", _vertexCandidateList.size());
- return;
+ int size = trackList.size();
+ for (int i = 0; i < size; ++i) {
+ ZvTrack iTrack = trackList.get(i);
+ for (int j = i + 1; j < size; ++j) {
+ ZvTrack 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();
}
/**
@@ -431,7 +392,6 @@
* 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<ZvVertex> fitVertices() {
@@ -464,4 +424,67 @@
}
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) {
+ ZvTrack iTrack = trackList.get(i);
+ double iTubeVal = iTrack.getTubeValue(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;
+ }
}