lcsim/src/org/lcsim/fit/twopointcircle
diff -u -r1.1 -r1.2
--- TwoPointCircleFit.java 29 Oct 2009 19:15:55 -0000 1.1
+++ TwoPointCircleFit.java 3 Nov 2009 23:52:49 -0000 1.2
@@ -12,6 +12,7 @@
private double _xc;
private double _yc;
private double _rc;
+ private boolean _cw;
private double _s1;
private double _s2;
@@ -21,13 +22,15 @@
* @param xc x coordinate of circle center
* @param yc y coordinate of circle center
* @param rc circle radius
+ * @param cw true if the closest hit to the DCA is in the clockwise direction
* @param s1 arc length from the DCA to the first hit
* @param s2 arc length from the DCA to the second hit
*/
- public TwoPointCircleFit(double xc, double yc, double rc, double s1, double s2) {
+ public TwoPointCircleFit(double xc, double yc, double rc, boolean cw, double s1, double s2) {
_xc = xc;
_yc = yc;
_rc = rc;
+ _cw = cw;
_s1 = s1;
_s2 = s2;
}
@@ -60,6 +63,16 @@
}
/**
+ * Return true if the closest hit to the DCA is in a clockwise direction
+ * relative to the the DCA
+ *
+ * @return true if hits are cw
+ */
+ public boolean cw() {
+ return _cw;
+ }
+
+ /**
* Return the arc length to the first hit.
*
* @return arc length to the first hit
lcsim/src/org/lcsim/fit/twopointcircle
diff -u -r1.1 -r1.2
--- TwoPointCircleFitter.java 29 Oct 2009 19:15:55 -0000 1.1
+++ TwoPointCircleFitter.java 3 Nov 2009 23:52:49 -0000 1.2
@@ -32,10 +32,11 @@
*/
public class TwoPointCircleFitter {
private double _rmin;
- private double _eps = 1e-6;
private List<TwoPointCircleFit> _circlefits;
private TwoPointLineFit _linefit;
private boolean _debug = true;
+ private static double _eps = 1e-6;
+ private static double twopi = 2. * Math.PI;
/**
* Constructor specifying a minimum radius. If the minimum radius is >0,
@@ -141,6 +142,9 @@
double r1dotr2 = x1*x2 + y1*y2;
double term1 = -ipinf * (r1dotr2 - dmaxsq) / denom;
double term2 = Math.abs(dmax * Math.sqrt((r1sq - dmaxsq) * (r2sq - dmaxsq)) / denom);
+
+ // Make sure alpha[0] is the solution with the largest circle radius
+ if (term1 < 0.) term2 = -term2;
alpha[0] = term1 + term2;
alpha[1] = term1 - term2;
}
@@ -157,11 +161,19 @@
if (i == 0 && !sltrack) return false;
// If we get here, either the first solution passed the pT cut or we are consistent with a
- // straight-line track. Find the circle with the minimum radius
+ // straight-line track. Find the circle with the minimum radius and the new alpha
rcurv = _rmin;
double newalpha = Math.sqrt(_rmin*_rmin - 0.25 * u*u);
- if (alpha[i] > 0) alpha[i] = newalpha;
- else alpha[i] = -newalpha;
+
+ // Figure out the sign of alpha. If we are consistent with a
+ // straight-line track, then the sign should be opposite for the
+ // second (low radius) solution. If we are not consistent with
+ // a straight-line track, the second solution should have the same
+ // sign as the first solution.
+ int isign = 1;
+ if (i == 1 && sltrack && alpha[0] >= 0.) isign = -1;
+ if (i == 1 && !sltrack && alpha[0] < 0.) isign = -1;
+ alpha[i] = isign * newalpha;
}
// Find the center of the circle
@@ -187,12 +199,37 @@
}
// Find the x-y arc lengths to hits 1 and 2
- double d1 = Math.sqrt((x1-x0)*(x1-x0) + (y1-y0)*(y1-y0));
- double s1 = Math.asin(0.5 * d1 / rcurv);
- double d2 = Math.sqrt((x2-x0)*(x2-x0) + (y2-y0)*(y2-y0));
- double s2 = Math.asin(0.5 * d2 / rcurv);
+ // First find azimuthal angles for the dca and hit positions relative to the circle center
+ double phi0 = Math.atan2(y0-yc, x0-xc);
+ double phi1 = Math.atan2(y1-yc, x1-xc);
+ double phi2 = Math.atan2(y2-yc, x2-xc);
+
+ // Find the angle between the hits and the DCA under the assumption that |dphi| < pi
+ double dphi1 = phi1 - phi0;
+ if (dphi1 > Math.PI) dphi1 -= twopi;
+ if (dphi1 < -Math.PI) dphi1 += twopi;
+ double dphi2 = phi2 - phi0;
+ if (dphi2 > Math.PI) dphi2 -= twopi;
+ if (dphi2 < -Math.PI) dphi2 += twopi;
+
+ // Find the hit closest to the DCA and use this hit to determine the circle "direction"
+ boolean cw;
+ if (Math.abs(dphi1) < Math.abs(dphi2)) cw = dphi1 < 0.;
+ else cw = dphi2 < 0.;
+
+ // Find the arc lengths to points 1 and 2
+ double s1 = -dphi1 * rcurv;
+ double s2 = -dphi2 * rcurv;
+ if (!cw) {
+ s1 = -s1;
+ s2 = -s2;
+ }
- _circlefits.add(new TwoPointCircleFit(xc, yc, rcurv, s1, s2));
+ // Fix the case when dphi1 & dphi2 have opposite signs, making an arc length negative
+ if (s1 < 0.) s1 += twopi * rcurv;
+ if (s2 < 0.) s2 += twopi * rcurv;
+
+ _circlefits.add(new TwoPointCircleFit(xc, yc, rcurv, cw, s1, s2));
}
// If we are consistent with a straight-line track, calculate the straight-line
@@ -202,9 +239,26 @@
// Find the distances from the straight-line DCA to the hits
double s1 = (x1*ux + y1*uy);
double s2 = (x2*ux + y2*uy);
- double phi = Math.atan2(y2-y1, x2-x1);
- _linefit = new TwoPointLineFit(phi, ipinf, s1, s2);
+ // If these distances have opposite sign, we violate causality for a hit originating at the DCA
+ if (s1*s2 < 0.) {
+
+ // Calculate the parameters of the line fit
+ double x0 = x1 - s1 * ux;
+ double y0 = y1 - s1 * uy;
+ double phi = Math.atan2(uy, ux);
+
+ // Flip the signs of the path lengths and direction to make path lengths positive
+ if (s1 < 0.) {
+ s1 = -s1;
+ s2 = -s2;
+ phi = phi + Math.PI;
+ if (phi > Math.PI) phi += twopi;
+ }
+
+ // Save a new TwoPointLineFit
+ _linefit = new TwoPointLineFit(x0, y0, phi, s1, s2);
+ }
// We should always find a valid circle above for this case - check to make sure
if (_debug & _circlefits.size() == 0) throw new RuntimeException("No circle found for hits consistent with infinite momentum");
lcsim/src/org/lcsim/fit/twopointcircle
diff -u -r1.1 -r1.2
--- TwoPointLineFit.java 29 Oct 2009 19:15:55 -0000 1.1
+++ TwoPointLineFit.java 3 Nov 2009 23:52:49 -0000 1.2
@@ -1,55 +1,64 @@
/*
* Encapsulate parameters for a straight line through two hits. The straight line
- * is parameterized as an azimuthal angle and an impact parameter, where the angle
- * defines the direction travelling from point 1 to point 2, the impact parameter
- * is positive if r1 x r2 is in the z direction, and s1 and s2 are the distances
- * to points 1 and 2 from the point of closest approach.
+ * is parameterized by specifying a point of closest approach and the azimuthal
+ * angle of the line.
*/
package org.lcsim.fit.twopointcircle;
/**
- * @author richp
+ * @author Richard Partridge
*/
public class TwoPointLineFit {
+ private double _x0;
+ private double _y0;
private double _phi0;
- private double _d0;
private double _s1;
private double _s2;
/**
* Fully qualified constructor.
*
+ * @param x0 x coordinate of point of closest approach
+ * @param y0 y coordinate of point of closest approach
* @param phi0 angle of line segment
- * @param d0 distance of closest approach
* @param s1 distance to point 1
* @param s2 distance to point 2
*/
- public TwoPointLineFit(double phi0, double d0, double s1, double s2) {
+ public TwoPointLineFit(double x0, double y0, double phi0, double s1, double s2) {
+ _x0 = x0;
+ _y0 = y0;
_phi0 = phi0;
- _d0 = d0;
_s1 = s1;
_s2 = s2;
}
/**
- * Return the azimuthal angle of the line segment going from point 1 to point 2
+ * Return the x coordinate of the point of closest approach
*
- * @return azimuthal angle
+ * @return x coordinate
*/
- public double phi0() {
- return _phi0;
+ public double x0() {
+ return _x0;
}
/**
- * Return the distance of closest approach for the line containing the two
- * points. The DCA is positive if r1 x r2 points in the +z direction.
+ * Return the y coordinate of the point of closest approach
*
- * @return distance of closest approach
+ * @return y coordinate
*/
- public double d0() {
- return _d0;
+ public double y0() {
+ return _y0;
+ }
+
+ /**
+ * Return the azimuthal angle of the line segment going from point 1 to point 2
+ *
+ * @return azimuthal angle
+ */
+ public double phi0() {
+ return _phi0;
}
/**
