 |
 |
| Commit in lcsim/src/org/lcsim/fit/helicaltrack on MAIN | |||
| HitUtils.java | +44 | -34 | 1.11 -> 1.12 |
Change the stereo angle calculation to check the quadrant of the dot-product.
diff -u -r1.11 -r1.12 --- HitUtils.java 16 Nov 2010 16:23:18 -0000 1.11 +++ HitUtils.java 22 May 2012 22:08:46 -0000 1.12 @@ -27,11 +27,12 @@
*/
public class HitUtils {
private static double _eps = 1.0e-6;
-
+ private boolean _debug = false; +
/** Creates a new instance of HitUtils */
public HitUtils() {
}
-
+
/**
* Turn a pixel hit into a pseudo-strip hit. Used for ZSegment hits that
* include a pixel hit. If the pixel is in an endcap disk, we need to
@@ -44,7 +45,7 @@
*/
public static HelicalTrack2DHit PixelToStrip(HelicalTrackHit hit, Map<HelicalTrackHit, Double> smap,
Map<HelicalTrackHit, MultipleScatter> msmap, HelicalTrackFit helix, double tolerance) {
-
+
// Take the strip to be span the interval +- tolerance * dz about the pixel hit position
double dz = zres(hit, msmap, helix);
double zmin = hit.z() - tolerance * dz;
@@ -56,14 +57,14 @@
hit.Detector(), hit.Layer(), hit.BarrelEndcapFlag(), zmin, zmax);
// Save the path length for the strip hit
smap.put(striphit, smap.get(hit));
-
+
return striphit;
}
-
+
public static Hep3Vector StripCenter(HelicalTrackStrip strip) {
return VecOp.add(strip.origin(), VecOp.mult(strip.umeas(), strip.u()));
}
-
+
public static SymmetricMatrix StripCov(HelicalTrackStrip strip) {
SymmetricMatrix cov = new SymmetricMatrix(3);
Hep3Vector pos = StripCenter(strip);
@@ -76,7 +77,7 @@
cov.setElement(1, 1, x*x * du*du / r2);
return cov;
}
-
+
public static SymmetricMatrix PixelCov(double x, double y, double drphi, double dz) {
SymmetricMatrix cov = new SymmetricMatrix(3);
double r2 = x*x + y*y;
@@ -86,7 +87,7 @@
cov.setElement(2, 2, dz);
return cov;
}
-
+
/**
* Find the effective z uncertainty to use in the s-z line fit. Include
* both resolution and multiple scattering contributions. Endcap disk
@@ -101,11 +102,11 @@
*/
public static double zres(HelicalTrackHit hit, Map<HelicalTrackHit, MultipleScatter> msmap,
HelicalTrackFit helix) {
-
+
// Get the multiple scattering uncertainty (if any)
double dz_ms = 0.;
if (msmap.containsKey(hit)) dz_ms = msmap.get(hit).dz();
-
+
// Barrels and disks are treated differently
if (hit.BarrelEndcapFlag() == BarrelEndcapFlag.BARREL) {
// For barrel hits, take the resolution uncertainty from the hit covariance matrix
@@ -120,32 +121,31 @@
if (helix != null) {
// Don't use the helix slope if the magnitude of the slope is smaller than its uncertainty
if (Math.abs(helix.slope()) > helix.getSlopeError()) slope = helix.slope();
-
+
}
// Take the resolution uncertainty to be dr * |slope|
double dzres = hit.dr() * Math.abs(slope);
// Combine resolution and multiple scattering uncertainties in quadrature
return Math.sqrt(dzres*dzres + dz_ms*dz_ms);
- }
+ }
}
-
+
/**
* Return the hit position assuming the track originated at the origin.
* @return hit position
*/
public static Hep3Vector PositionFromOrigin(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
-
+
// Assume the particle is coming from the origin, so x2 = gamma * x1
// gamma = Origin2 . w_hat / Origin1 . w_hat
double gamma = VecOp.dot(strip2.origin(), strip2.w()) / NonZeroDotProduct(strip1.origin(), strip1.w());
-
+
// Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
- double salpha = VecOp.dot(strip1.v(), strip2.u()); -
+ double salpha = getSinAlpha(strip1, strip2);
// Calculate the hit locations for v = 0: p = Origin + u * u_hat
Hep3Vector p1 = StripCenter(strip1);
Hep3Vector p2 = StripCenter(strip2);
-
+
// dp = (p2 - gamma * p1)
Hep3Vector dp = VecOp.sub(p2, VecOp.mult(gamma, p1));
@@ -157,9 +157,9 @@
// Position of hit on strip 1: r1 = p1 + v1 * v1_hat
Hep3Vector r1 = VecOp.add(p1, VecOp.mult(v1, strip1.v()));
// Take position to be the midpoint of r1 and r2: r = 0.5 * (1 + gamma) * r1
- return VecOp.mult(0.5 * (1 + gamma), r1);
+ return VecOp.mult(0.5 * (1 + gamma), r1);
}
-
+
/**
* Return the covariance matrix assuming the track originated from the
* origin.
@@ -172,7 +172,7 @@
// Calculate the seperation between the sensor planes
double separation = SensorSeperation(strip1, strip2);
// Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
- double salpha = v1Dotu2(strip1, strip2);
+ double salpha = getSinAlpha(strip1, strip2);
// Calculate the scale factor: factor = (1 + gamma)^2 / 4 * sin^2(alpha)
double factor = (1 + gamma)*(1 + gamma) / (4. * salpha*salpha);
// Calculate v * v^T for strips 1 and 2
@@ -189,7 +189,7 @@
// Don't let dv by greater than the strip length / sqrt(12)
double dv1 = Math.min(dv, (strip1.vmax()-strip1.vmin()) / Math.sqrt(12.));
double dv2 = Math.min(dv, (strip2.vmax()-strip2.vmin()) / Math.sqrt(12.));
- // Calculate the covariance matrix.
+ // Calculate the covariance matrix.
// From resolution: cov = factor * (v2 * v2^T * du1^2 + v1 * v1^T * du2^2)
// From direction: + (v1 * v1^T * (dv1/2)^2 + v2 * v2^T * (dv2/2)^2
Matrix cov1 = MatrixOp.mult(factor * du2*du2 + 0.25 * dv1*dv1, v1v1t);
@@ -197,7 +197,7 @@
Matrix cov = MatrixOp.add(cov1, cov2);
return new SymmetricMatrix(cov);
}
-
+
/**
* Return the hit position given the track direction.
* @param trkdir TrackDirection object containing direction and derivatives
@@ -207,7 +207,7 @@
HelicalTrackStrip strip2) {
// Get the track direction unit vector
Hep3Vector dir = trkdir.Direction();
- double salpha = v1Dotu2(strip1, strip2);
+ double salpha = getSinAlpha(strip1, strip2);
// Gamma is the distance the particle travels between the two sensors: gamma = separation / (what . dir)
double gamma = SensorSeperation(strip1, strip2) / NonZeroDotProduct(strip1.w(), dir);
Hep3Vector p1 = StripCenter(strip1);
@@ -221,7 +221,7 @@
// Take position to be the midpoint of x1 and x2: r = r1 + 0.5 * gamma * dir
return VecOp.add(r1, VecOp.mult(0.5 * gamma, dir));
}
-
+
/**
* Return the covariance matrix given a track direction and helix
@@ -242,7 +242,7 @@
// Calculate phat . w
double pdotw = NonZeroDotProduct(dir, strip1.w());
// salpha is the sin of the stereo angle
- double salpha = v1Dotu2(strip1, strip2);
+ double salpha = getSinAlpha(strip1, strip2);
// Calculate the scale factor: _separation / 2 * sin(alpha) * (phat . w)^2
double factor = SensorSeperation(strip1, strip2) / (2 * salpha * pdotw*pdotw);
// Create the matrix of position derivatives: d_i,j = dr_i / dphat_j
@@ -250,7 +250,7 @@
Matrix v1 = v2m(strip1.v());
Matrix v2 = v2m(strip2.v());
Matrix d = MatrixOp.mult(factor, MatrixOp.add(MatrixOp.mult(v1, MatrixOp.transposed(pcrossv2)),
- MatrixOp.mult(v2, MatrixOp.transposed(pcrossv1))));
+ MatrixOp.mult(v2, MatrixOp.transposed(pcrossv1))));
Matrix dh = MatrixOp.mult(d, dirderiv);
// Construct the transpose of dh
Matrix dht = MatrixOp.transposed(dh);
@@ -267,14 +267,14 @@
// Convert to a symmetric matrix
return new SymmetricMatrix(cov);
}
-
+
public static double UnmeasuredCoordinate(TrackDirection trkdir, HelicalTrackStrip strip1,
HelicalTrackStrip strip2) {
// Get the track direction unit vector
Hep3Vector dir = trkdir.Direction();
// Gamma is the distance the particle travels between the two sensors: gamma = separation / (what . dir)
double gamma = SensorSeperation(strip1, strip2) / NonZeroDotProduct(strip1.w(), dir);
- double salpha = v1Dotu2(strip1, strip2);
+ double salpha = getSinAlpha(strip1, strip2);
Hep3Vector p1 = StripCenter(strip1);
Hep3Vector p2 = StripCenter(strip2);
// dp = p2 - (p1 + gamma * dir)
@@ -282,7 +282,7 @@
// Unmeasured coordinate v1: v1 = (dp . u2hat) / sin(alpha)
return VecOp.dot(dp, strip2.u()) / salpha;
}
-
+
/**
* Calculate the uncertainty in the unmeasured coordinate v1.
* @param trkdir track direction
@@ -300,7 +300,7 @@
double pdotw = NonZeroDotProduct(dir, strip1.w());
// Calculate phat x v2
Hep3Vector pcrossv2 = VecOp.cross(dir, strip2.v());
- double salpha = v1Dotu2(strip1, strip2);
+ double salpha = getSinAlpha(strip1, strip2);
// Calculate the scale factor: separation / sin(alpha) * (phat . w)^2
double factor = SensorSeperation(strip1, strip2) / (salpha * pdotw*pdotw);
// The matrix d^T is a row vector given by factor * (phat x v2)
@@ -318,17 +318,27 @@
double dvsq = (Math.pow(u1dotu2 * strip1.du(), 2) + Math.pow(strip2.du(), 2))/ (salpha*salpha);
return Math.sqrt(dvsq + cov.e(0, 0));
}
-
+
public static double SensorSeperation(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
// Calculate the seperation between the sensor planes
return VecOp.dot(strip1.w(), VecOp.sub(strip2.origin(), strip1.origin()));
}
-
+
public static double v1Dotu2(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
// Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
return VecOp.dot(strip1.v(), strip2.u());
}
-
+
+ public static double getSinAlpha(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
+ // Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
+ double salpha = v1Dotu2(strip1, strip2);
+ //get cos(alpha) and check if the meaurement directions are ~parallel or flipped
+ double calpha = VecOp.dot(strip1.u(), strip2.u());
+ if (calpha < 0)
+ salpha = -salpha;
+ return salpha;
+ }
+
private static double NonZeroDotProduct(Hep3Vector v1, Hep3Vector v2) {
double cth = VecOp.dot(v1, v2);
if (Math.abs(cth) < _eps) {
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1