lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.6 -r1.7
--- HelixUtils.java 17 Jun 2008 00:40:49 -0000 1.6
+++ HelixUtils.java 13 Oct 2008 01:05:58 -0000 1.7
@@ -30,25 +30,7 @@
* @version 1.0
*/
public class HelixUtils {
- private CircleFit _cfit;
- private HelicalTrackFit _helix;
- private double _dca;
- private double _phi0;
- private double _omega;
- private double _RC;
- private double _z0;
- private double _slope;
- private double _xr;
- private double _yr;
- private double _xc;
- private double _yc;
- private double _x0;
- private double _y0;
- private double _sphi0;
- private double _cphi0;
- private double _sth;
- private double _cth;
- private double _minslope = 1.e-6;
+ private static double _minslope = 1.e-6;
/**
* Creates a new instance of HelixUtils.
@@ -63,11 +45,9 @@
* @param hit2 second hit
* @return path length from hit1 to hit2
*/
- public double PathLength(CircleFit cfit, HelicalTrackHit hit1, HelicalTrackHit hit2) {
- // Cache the circle parameters
- CircleCache(cfit);
+ public static double PathLength(CircleFit cfit, HelicalTrackHit hit1, HelicalTrackHit hit2) {
// Return the path length between hit1 and hit2
- return PathCalc(hit1.x(), hit1.y(), hit2.x(), hit2.y());
+ return PathCalc(xc(cfit), yc(cfit), RC(cfit), hit1.x(), hit1.y(), hit2.x(), hit2.y());
}
/**
@@ -76,11 +56,9 @@
* @param hit hit to be used for the path length calculation
* @return path length from the DCA to the hit
*/
- public double PathLength(CircleFit cfit, HelicalTrackHit hit) {
- // Cache the circle parameters
- CircleCache(cfit);
+ public static double PathLength(CircleFit cfit, HelicalTrackHit hit) {
// Return the path length from the DCA
- return PathCalc(_x0, _y0, hit.x(), hit.y());
+ return PathCalc(xc(cfit), yc(cfit), RC(cfit), x0(cfit), y0(cfit), hit.x(), hit.y());
}
/**
@@ -89,11 +67,9 @@
* @param hit hit to be used for the path length calculation
* @return path length from the DCA to the hit
*/
- public double PathLength(HelicalTrackFit helix, HelicalTrackHit hit) {
- // Cache the helix parameters
- HelixCache(helix);
+ public static double PathLength(HelicalTrackFit helix, HelicalTrackHit hit) {
// Return the path length from the DCA
- return PathCalc(_x0, _y0, hit.x(), hit.y());
+ return PathCalc(helix.xc(), helix.yc(), helix.R(), helix.x0(), helix.y0(), hit.x(), hit.y());
}
/**
@@ -102,13 +78,11 @@
* @param z location of z-plane
* @return path length from the DCA to the z-plane
*/
- public double PathToZPlane(HelicalTrackFit helix, double z) {
- // Cache the helix parameters
- HelixCache(helix);
+ public static double PathToZPlane(HelicalTrackFit helix, double z) {
// Find the z distace between the DCA and the z plane and calculate the x-y path length
- double zdist = z - _z0;
- double safeslope = _slope;
- if (Math.abs(_slope) < _minslope) safeslope = _minslope * Math.signum(_slope);
+ double zdist = z - helix.z0();
+ double safeslope = helix.slope();
+ if (Math.abs(safeslope) < _minslope) safeslope = _minslope * Math.signum(safeslope);
return zdist / safeslope;
}
@@ -120,28 +94,29 @@
* @param mxint Maximum number of intersections
* @return list of path lengths
*/
- public List<Double> PathToCylinder(HelicalTrackFit helix, double r, double smax, int mxint) {
- // Cache the helix parameters
- HelixCache(helix);
+ public static List<Double> PathToCylinder(HelicalTrackFit helix, double r, double smax, int mxint) {
// Create a list to hold the path lengths
List<Double> pathlist = new ArrayList<Double>();
+ // Retrieve helix dca and RC
+ double dca = helix.dca();
+ double RC = helix.R();
// Find cos(dphi) for dphi from the point of closest approach to the first crossing
- double cdphi = 1 - (r * r - _dca * _dca) / (2. * _RC * (_RC - _dca));
+ double cdphi = 1 - (r * r - dca * dca) / (2. * RC * (RC - dca));
// Make sure we have a physical intersection
if (Math.abs(cdphi) < 1.) {
// Calculate dphi and the path length to the first crossing
double dphi = Math.acos(cdphi);
- Double s = dphi * Math.abs(_RC);
+ Double s = dphi * Math.abs(RC);
// Loop over crossings until we exceed one of the limits
while (s < smax && pathlist.size() < mxint) {
// Add this odd-numbered crossing to the list
pathlist.add(s);
// Advance to the next even-numbered crossing
- s += 2. * (Math.PI - dphi) * Math.abs(_RC);
+ s += 2. * (Math.PI - dphi) * Math.abs(RC);
// Check to see if we should add it
if (s < smax && pathlist.size() < mxint) pathlist.add(s);
// Add this even-numbered crossing to the list
- s += 2. * dphi * Math.abs(_RC);
+ s += 2. * dphi * Math.abs(RC);
}
}
return pathlist;
@@ -154,15 +129,14 @@
* @param s path length to the desired point on helix
* @return direction unit vector
*/
- public Hep3Vector Direction(HelicalTrackFit helix, double s) {
- // Cache the helix parameters
- HelixCache(helix);
+ public static Hep3Vector Direction(HelicalTrackFit helix, double s) {
// Calculate the azimuthal direction
- double phi = _phi0 - s / _RC;
+ double phi = helix.phi0() - s / helix.R();
+ double sth = helix.sth();
// Calculate the components of the direction unit vector
- double ux = Math.cos(phi) * _sth;
- double uy = Math.sin(phi) * _sth;
- double uz = _cth;
+ double ux = Math.cos(phi) * helix.sth();
+ double uy = Math.sin(phi) * helix.sth();
+ double uz = helix.cth();
// Return the direction unit vector
return new BasicHep3Vector(ux, uy, uz);
}
@@ -176,7 +150,7 @@
* @param s path length specifying location on helix
* @return TrackDirection object for this point on the helix
*/
- public TrackDirection CalculateTrackDirection(HelicalTrackFit helix, double s) {
+ public static TrackDirection CalculateTrackDirection(HelicalTrackFit helix, double s) {
Hep3Vector dir = Direction(helix, s);
Matrix deriv = DirectionDerivates(helix, dir, s);
return new TrackDirection(dir, deriv);
@@ -188,29 +162,28 @@
* @param s path length
* @return point in space corresponding to the given path length
*/
- public Hep3Vector PointOnHelix(HelicalTrackFit helix, double s) {
- // Cache the helix parameters
- HelixCache(helix);
+ public static Hep3Vector PointOnHelix(HelicalTrackFit helix, double s) {
// Find the azimuthal direction at this path length
- double phi = _phi0 - s / _RC;
+ double RC = helix.R();
+ double phi = helix.phi0() - s / RC;
// Calculate the position on the helix at this path length
- double x = _xc - _RC * Math.sin(phi);
- double y = _yc + _RC * Math.cos(phi);
- double z = _z0 + s * _slope;
+ double x = helix.xc() - RC * Math.sin(phi);
+ double y = helix.yc() + RC * Math.cos(phi);
+ double z = helix.z0() + s * helix.slope();
// Return the position as a Hep3Vector
return new BasicHep3Vector(x, y, z);
}
- private double PathCalc(double x1, double y1, double x2, double y2) {
+ private static double PathCalc(double xc, double yc, double RC, double x1, double y1, double x2, double y2) {
// Find the angle between these points measured wrt the circle center
- double phi1 = Math.atan2(y1 - _yc, x1 - _xc);
- double phi2 = Math.atan2(y2 - _yc, x2 - _xc);
+ double phi1 = Math.atan2(y1 - yc, x1 - xc);
+ double phi2 = Math.atan2(y2 - yc, x2 - xc);
double dphi = phi2 - phi1;
// Make sure dphi is in the valid range (-pi, pi)
if (dphi > Math.PI) dphi -= 2. * Math.PI;
if (dphi < -Math.PI) dphi += 2. * Math.PI;
// Return the arc length
- return -_RC * dphi;
+ return -RC * dphi;
}
/**
@@ -222,67 +195,47 @@
* @param s path length to the desired point on the helix
* @return direction derivatives
*/
- private Matrix DirectionDerivates(HelicalTrackFit helix, Hep3Vector u, double s) {
- // Cache the helix parameters
- HelixCache(helix);
+ private static Matrix DirectionDerivates(HelicalTrackFit helix, Hep3Vector u, double s) {
// Create the matrix that will hold the derivatives
BasicMatrix deriv = new BasicMatrix(3,5);
+ // Retrieve some helix info
+ double cphi0 = Math.cos(helix.phi0());
+ double sphi0 = Math.sin(helix.phi0());
+ double sth = helix.sth();
+ double cth = helix.cth();
+ double dca = helix.dca();
+ double omega = helix.curvature();
// Calculate the non-zero derivatives of the direction with respect to the helix parameters
- deriv.setElement(0, helix.curvatureIndex, (u.x() - _cphi0 * _sth) / _omega); // du_x / domega
- deriv.setElement(1, helix.curvatureIndex, (u.y() - _sphi0 * _sth) / _omega); // du_y / domega
- deriv.setElement(0, helix.dcaIndex, -_omega * _cphi0 * _sth); // du_x / ddca
- deriv.setElement(1, helix.dcaIndex, -_omega * _sphi0 * _sth); // du_y / ddca
- deriv.setElement(0, helix.phi0Index, -(1 - _dca * _omega) * _sphi0 * _sth); // du_x / dphi0
- deriv.setElement(1, helix.phi0Index, (1 - _dca * _omega) * _cphi0 * _sth); // du_y / dphi0
- deriv.setElement(0, helix.slopeIndex, -u.x() * _sth * _cth); // du_x / dslope
- deriv.setElement(1, helix.slopeIndex, -u.y() * _sth * _cth); // du_y / dslope
- deriv.setElement(2, helix.slopeIndex, _sth*_sth*_sth); // du_z / dslope
+ deriv.setElement(0, helix.curvatureIndex, (u.x() - cphi0 * sth) / omega); // du_x / domega
+ deriv.setElement(1, helix.curvatureIndex, (u.y() - sphi0 * sth) / omega); // du_y / domega
+ deriv.setElement(0, helix.dcaIndex, -omega * cphi0 * sth); // du_x / ddca
+ deriv.setElement(1, helix.dcaIndex, -omega * sphi0 * sth); // du_y / ddca
+ deriv.setElement(0, helix.phi0Index, -(1 - dca * omega) * sphi0 * sth); // du_x / dphi0
+ deriv.setElement(1, helix.phi0Index, (1 - dca * omega) * cphi0 * sth); // du_y / dphi0
+ deriv.setElement(0, helix.slopeIndex, -u.x() * sth * cth); // du_x / dslope
+ deriv.setElement(1, helix.slopeIndex, -u.y() * sth * cth); // du_y / dslope
+ deriv.setElement(2, helix.slopeIndex, sth*sth*sth); // du_z / dslope
// Return the derivative matrix
return deriv;
}
+
+ private static double RC(CircleFit cfit) {
+ return 1.0 / cfit.curvature();
+ }
- private void CircleCache(CircleFit cfit) {
- // Cache some calculations for this circle if we haven't already saved them
- if (_cfit != cfit) {
- _cfit = cfit;
- _helix = null;
- _dca = -_cfit.dca();
- _phi0 = _cfit.phi();
- _RC = 1.0 / _cfit.curvature();
- _xr = _cfit.xref();
- _yr = _cfit.yref();
- _sphi0 = Math.sin(_phi0);
- _cphi0 = Math.cos(_phi0);
- _xc = _xr + (_RC - _dca) * _sphi0;
- _yc = _yr - (_RC - _dca) * _cphi0;
- _x0 = _xr - _dca * _sphi0;
- _y0 = _yr + _dca * _cphi0;
- }
- return;
+ private static double xc(CircleFit cfit) {
+ return cfit.xref() + (RC(cfit) - cfit.dca()) * Math.sin(cfit.phi());
}
- private void HelixCache(HelicalTrackFit helix) {
- // Cache some calculations for this helix if we haven't already saved them
- if (_helix != helix) {
- _cfit = null;
- _helix = helix;
- _dca = _helix.dca();
- _phi0 = _helix.phi0();
- _omega = _helix.curvature();
- _RC = 1.0 / _omega;
- _z0 = _helix.z0();
- _slope = _helix.slope();
- _xr = 0.;
- _yr = 0.;
- _sphi0 = Math.sin(_phi0);
- _cphi0 = Math.cos(_phi0);
- _xc = _xr + (_RC - _dca) * _sphi0;
- _yc = _yr - (_RC - _dca) * _cphi0;
- _x0 = _xr - _dca * _sphi0;
- _y0 = _yr + _dca * _cphi0;
- _sth = _helix.sth();
- _cth = _helix.cth();
- }
- return;
+ private static double yc(CircleFit cfit) {
+ return cfit.yref() - (RC(cfit) - cfit.dca()) * Math.cos(cfit.phi());
+ }
+
+ private static double x0(CircleFit cfit) {
+ return cfit.xref() - cfit.dca() * Math.sin(cfit.phi());
+ }
+
+ private static double y0(CircleFit cfit) {
+ return cfit.yref() + cfit.dca() * Math.cos(cfit.phi());
}
}
\ No newline at end of file