lcsim/src/org/lcsim/recon/tracking/seedtracker
diff -u -r1.10 -r1.11
--- FastCheck.java 4 Sep 2009 23:17:01 -0000 1.10
+++ FastCheck.java 6 Nov 2009 00:29:40 -0000 1.11
@@ -2,14 +2,18 @@
import org.lcsim.constants.Constants;
import org.lcsim.event.MCParticle;
+import org.lcsim.fit.threepointcircle.CircleFit;
import org.lcsim.fit.helicaltrack.HelicalTrack2DHit;
+import org.lcsim.fit.helicaltrack.HelicalTrack3DHit;
import org.lcsim.fit.helicaltrack.HelicalTrackCross;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
-import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
-import org.lcsim.fit.helicaltrack.HelixParamCalculator;
import org.lcsim.fit.helicaltrack.HelixUtils;
import org.lcsim.fit.helicaltrack.TrackDirection;
+import org.lcsim.fit.threepointcircle.ThreePointCircleFitter;
+import org.lcsim.fit.twopointcircle.TwoPointCircleFit;
+import org.lcsim.fit.twopointcircle.TwoPointCircleFitter;
+import org.lcsim.fit.twopointcircle.TwoPointLineFit;
import org.lcsim.recon.tracking.seedtracker.diagnostic.ISeedTrackerDiagnostics;
/**
@@ -18,12 +22,15 @@
*/
public class FastCheck {
-// private HelicalTrackHit _hit1cache;
private SeedStrategy _strategy;
private double _bfield;
private double _RMin;
private double _dMax;
private double _z0Max;
+ private double _nsig;
+ private TwoPointCircleFitter _cfit2;
+ private ThreePointCircleFitter _cfit3;
+ private static double twopi = 2. * Math.PI;
private double _eps = 1.0e-6;
private ISeedTrackerDiagnostics _diag;
@@ -36,114 +43,20 @@
_RMin = strategy.getMinPT() / (Constants.fieldConversion * bfield);
_dMax = strategy.getMaxDCA();
_z0Max = strategy.getMaxZ0();
- }
-
- public boolean CheckHitPair(HelicalTrackHit hit1, HelicalTrackHit hit2, SeedCandidate seed) {
-
- // Check if seed and hits are all from a single "True Seed"
- boolean match = false;
- if (_diag != null && seed.isTrueSeed()) {
- for (MCParticle mcp : seed.getMCParticles()) {
- if (hit1.getMCParticles().contains(mcp) && hit2.getMCParticles().contains(mcp)) {
- match = true;
- break;
- }
- }
- }
+ _nsig = Math.sqrt(strategy.getMaxChisq());
- // Correct the hit positions for stereo hits
- CorrectHitPosition(hit1, seed);
- CorrectHitPosition(hit2, seed);
+ // Instantiate the two point circle fitter for this minimum radius, maximum DCA
+ _cfit2 = new TwoPointCircleFitter(_RMin);
- // Get the polar coordinates for the hits
- double r1 = hit1.r();
- double r2 = hit2.r();
- double phi1 = hit1.phi();
- double phi2 = hit2.phi();
-
- // Get the maximum deviation in phi between these two radii
- double dphimx = dphimax(r1, r2);
-
- // Calculate the difference in azimuthal angle between these two hits, which should always be positive
- double dphi = phidif(phi2, phi1);
-
- // Check that the difference in azimuthal angle is less than the difference in maximum phi deviations
- boolean phiOK = dphi <= dphimx;
-
- if (!phiOK) {
- if (match) {
- System.out.println("Bad Phi check");
-// System.out.println("dphi: " + dphi + " dphimx: " + dphimx);
-// System.out.println("r1: " + r1 + " r2: " + r2 + " phi1: " + phi1 + " phi2: " + phi2);
- }
- return false;
- }
-
- // Get the z limits for the hits
- double zlen1 = 0.;
- double zlen2 = 0.;
- if (hit1 instanceof HelicalTrack2DHit) {
- zlen1 = ((HelicalTrack2DHit) hit1).zlen();
- }
- if (hit2 instanceof HelicalTrack2DHit) {
- zlen2 = ((HelicalTrack2DHit) hit2).zlen();
- }
- double zmin1 = hit1.z() - 0.5 * zlen1;
- double zmax1 = zmin1 + zlen1;
- double zmin2 = hit2.z() - 0.5 * zlen2;
- double zmax2 = zmin2 + zlen2;
-
- // Check the z0 limits using the minimum path lengths
- boolean zOK = checkz0(smin(r1), smax(r1), zmin1, zmax1, smin(r2), smax(r2), zmin2, zmax2);
-
- if (!zOK) {
- if (match) {
- for (MCParticle mcp : seed.getMCParticles()) {
-// HelixParamCalculator helmc = new HelixParamCalculator(mcp, seed.getBField());
-// double z1 = 0.5 * (zmin1 + zmax1);
-// double z2 = 0.5 * (zmin2 + zmax2);
-// double s1true = (z1 - helmc.getZ0()) / helmc.getSlopeSZPlane();
-// double s2true = (z2 - helmc.getZ0()) / helmc.getSlopeSZPlane();
-// double arg1 = s1true * helmc.getMCOmega() / 2.;
-// double x1true = helmc.getX0() + s1true * Math.sin(arg1) * Math.cos(helmc.getPhi0()-arg1) / arg1;
-// double y1true = helmc.getY0() + s1true * Math.sin(arg1) * Math.sin(helmc.getPhi0()-arg1) / arg1;
-// double arg2 = s2true * helmc.getMCOmega() / 2.;
-// double x2true = helmc.getX0() + s2true * Math.sin(arg2) * Math.cos(helmc.getPhi0()-arg2) / arg2;
-// double y2true = helmc.getY0() + s2true * Math.sin(arg2) * Math.sin(helmc.getPhi0()-arg2) / arg2;
- // System.out.println("True s1: "+s1true+" x1: "+x1true+" y1: "+y1true);
- // System.out.println("True s2: "+s2true+" x2: "+x2true+" y2: "+y2true);
- }
- System.out.println("******* Bad Z Check *******");
-// System.out.println(" x1: "+hit1.getCorrectedPosition().x()+
-// " y1: "+hit1.getCorrectedPosition().y()+
-// " x2: "+hit2.getCorrectedPosition().x()+
-// " y2: "+hit2.getCorrectedPosition().y());
-// System.out.println("s1min: " + smin(r1) + " s1max: " + smax(r1) + " zmin1: " + zmin1 +
-// " zmax1: " + zmax1);
-// System.out.println("s2min: " + smin(r2) + " s2max: " + smax(r2) + " zmin2: " + zmin2 +
-// " zmax2: " + zmax2);
-// System.out.println("Hit1: "+hit1.toString());
-// System.out.println("Hit 2: "+hit2.toString());
-// if (hit2 instanceof HelicalTrackCross) {
-// HelicalTrackCross c2 = (HelicalTrackCross) hit2;
-// for (HelicalTrackStrip strip : c2.getStrips()) {
-// System.out.println("u: "+strip.u().toString());
-// System.out.println("v: "+strip.v().toString());
-// System.out.println("origin: "+strip.origin().toString());
-// System.out.println("umeas: "+strip.umeas()+" vmin: "+strip.vmin()+" vmax: "+strip.vmax());
-// }
-
-// }
- }
- }
- return zOK;
+ // Instantiate the three point circle fitter
+ _cfit3 = new ThreePointCircleFitter();
}
public boolean CheckHitSeed(HelicalTrackHit hit, SeedCandidate seed) {
// Check the hit against each hit in the seed
for (HelicalTrackHit hit2 : seed.getHits()) {
- if (!CheckHitPair(hit, hit2, seed)) return false;
+ if (!TwoPointCircleCheck(hit, hit2, seed)) return false;
}
return true;
}
@@ -216,30 +129,6 @@
if (match) {
for (MCParticle mcp : seed.getMCParticles()) {
System.out.println("z sector error - p: "+mcp.getMomentum().toString());
-// System.out.println("Sector ID: "+sector.Identifier()+" z0max: "+_z0Max);
-// System.out.println("zlo: "+sector.zlo()+"zhi: "+sector.zhi());
-// System.out.println("smin1: "+smin1+" smax1: "+smax1+" zmin1: "+zmin+" zmax1: "+zmax);
-// System.out.println("smin2: "+smin2+" smax2: "+smax2+" zmin2: "+zmin2+" zmax2: "+zmax2);
-// System.out.println("Hit: "+hit.toString());
-// HelixParamCalculator helmc = new HelixParamCalculator(mcp, seed.getBField());
-// System.out.println("Helix z: "+helmc.getZ0()+" slope: "+helmc.getSlopeSZPlane());
-// double z1 = 0.5 * (zmin + zmax);
-// double z2 = 0.5 * (zmin2 + zmax2);
-// double s1true = (z1 - helmc.getZ0()) / helmc.getSlopeSZPlane();
-// double s2true = (z2 - helmc.getZ0()) / helmc.getSlopeSZPlane();
-// double arg1 = s1true * helmc.getMCOmega() / 2.;
-// double x1true = helmc.getX0() + s1true * Math.sin(arg1) * Math.cos(helmc.getPhi0()-arg1) / arg1;
-// double y1true = helmc.getY0() + s1true * Math.sin(arg1) * Math.sin(helmc.getPhi0()-arg1) / arg1;
-// double arg2 = s2true * helmc.getMCOmega() / 2.;
-// double x2true = helmc.getX0() + s2true * Math.sin(arg2) * Math.cos(helmc.getPhi0()-arg2) / arg2;
-// double y2true = helmc.getY0() + s2true * Math.sin(arg2) * Math.sin(helmc.getPhi0()-arg2) / arg2;
-// System.out.println("True s1: "+s1true+" x1: "+x1true+" y1: "+y1true);
-// System.out.println("True s2: "+s2true+" x2: "+x2true+" y2: "+y2true);
-// for (HelicalTrackHit secthit : sector.Hits()) {
-// if (secthit.getMCParticles().contains(mcp)) {
-// System.out.println("Found true hit in sector: "+secthit.toString());
-// }
-// }
}
}
return false;
@@ -288,6 +177,206 @@
return zOK;
}
+ public boolean TwoPointCircleCheck(HelicalTrackHit hit1, HelicalTrackHit hit2, SeedCandidate seed) {
+
+ // Initialize the hit coordinates for an unknown track direction
+ CorrectHitPosition(hit1, seed);
+ CorrectHitPosition(hit2, seed);
+
+ // Check that hits are outside the maximum DCA
+ if (hit1.r() < _dMax || hit2.r() < _dMax) return false;
+
+ // Try to find a circle passing through the 2 hits and the maximum DCA
+ boolean success = _cfit2.FitCircle(hit1, hit2, _dMax);
+
+ // Check for success
+ if (!success) return false;
+
+ // Initialize the minimum/maximum arc lengths
+ double s1min = 1.0e99;
+ double s1max = -1.0e99;
+ double s2min = 1.0e99;
+ double s2max = -1.0e99;
+
+ // Loop over the circle fits and find the min/max arc lengths
+ for (TwoPointCircleFit fit : _cfit2.getCircleFits()) {
+ double s1 = fit.s1();
+ double s2 = fit.s2();
+ if (s1 < s1min) s1min = s1;
+ if (s1 > s1max) s1max = s1;
+ if (s2 < s2min) s2min = s2;
+ if (s2 > s2max) s2max = s2;
+ }
+
+ // If we are consistent with a straight-line fit, update the minimum s1, s2
+ TwoPointLineFit lfit = _cfit2.getLineFit();
+ if (lfit != null) {
+
+ // Find the distance from the DCA to the maximum DCA circle
+ double x0 = lfit.x0();
+ double y0 = lfit.y0();
+ double s0 = 0.;
+ double s0sq = _dMax*_dMax - (x0*x0 + y0*y0);
+ if (s0sq > _eps*_eps) s0 = Math.sqrt(s0sq);
+
+ // Update the minimum arc length to the distance from the DCA to the hit
+ s1min = lfit.s1() - s0;
+ s2min = lfit.s2() - s0;
+ }
+
+ // Now check for consistent hits in the s-z plane
+ // First get the z limits for the hits
+ double z1len = 0.;
+ double z2len = 0.;
+ if (hit1 instanceof HelicalTrack2DHit) {
+ z1len = ((HelicalTrack2DHit) hit1).zlen();
+ }
+ if (hit2 instanceof HelicalTrack2DHit) {
+ z2len = ((HelicalTrack2DHit) hit2).zlen();
+ }
+ double z1min = hit1.z() - 0.5 * z1len;
+ double z1max = z1min + z1len;
+ double z2min = hit2.z() - 0.5 * z2len;
+ double z2max = z2min + z2len;
+
+ // Check the z0 limits using the min/max path lengths
+ boolean zOK = checkz0(s1min, s1max, z1min, z1max, s2min, s2max, z2min, z2max);
+
+ // Done!
+ return zOK;
+ }
+
+ public boolean ThreePointHelixCheck(HelicalTrackHit hit1, HelicalTrackHit hit2, HelicalTrackHit hit3) {
+
+ // Setup for a 3 point circle fit
+ double p[][] = new double[3][2];
+ double[] pos;
+ double z[] = new double[3];
+ double dztot = 0.;
+ boolean zfirst = true;
+
+ // While not terribly elegant, code for speed
+ // Use calls that give uncorrected position and error
+ // Get the relevant variables for hit 1
+ pos = hit1.getPosition();
+ p[0][0] = pos[0];
+ p[0][1] = pos[1];
+ z[0] = pos[2];
+ if (hit1 instanceof HelicalTrack3DHit) dztot += _nsig * ((HelicalTrack3DHit) hit1).dz();
+ else {
+ zfirst = false;
+ if (hit1 instanceof HelicalTrack2DHit) dztot += ((HelicalTrack2DHit) hit1).zlen() / 2.;
+ else dztot += _nsig * Math.sqrt(hit1.getCovMatrix()[5]);
+ }
+
+ // Get the relevant variables for hit 2
+ pos = hit2.getPosition();
+ p[1][0] = pos[0];
+ p[1][1] = pos[1];
+ z[1] = pos[2];
+ if (hit2 instanceof HelicalTrack3DHit) dztot += _nsig * ((HelicalTrack3DHit) hit2).dz();
+ else {
+ zfirst = false;
+ if (hit2 instanceof HelicalTrack2DHit) dztot += ((HelicalTrack2DHit) hit2).zlen() / 2.;
+ else dztot += _nsig * Math.sqrt(hit2.getCovMatrix()[5]);
+ }
+
+ // Get the relevant variables for hit 3
+ pos = hit3.getPosition();
+ p[2][0] = pos[0];
+ p[2][1] = pos[1];
+ z[2] = pos[2];
+ if (hit3 instanceof HelicalTrack3DHit) dztot += _nsig * ((HelicalTrack3DHit) hit3).dz();
+ else {
+ zfirst = false;
+ if (hit3 instanceof HelicalTrack2DHit) dztot += ((HelicalTrack2DHit) hit3).zlen() / 2.;
+ else dztot += _nsig * Math.sqrt(hit3.getCovMatrix()[5]);
+ }
+
+ // Add multiple scattering error here - for now, just set it to 1 mm
+ dztot += 1.;
+
+ // Unless the three hits are all pixel hits, do the circle checks first
+ if (!zfirst) {
+ if (!TwoPointCircleCheck(hit1, hit3, null)) return false;
+ if (!TwoPointCircleCheck(hit2, hit3, null)) return false;
+ }
+
+ // Do the 3 point circle fit and check for success
+ boolean success = _cfit3.fit(p[0], p[1], p[2]);
+ if (!success) return false;
+
+ // Retrieve the circle parameters
+ CircleFit circle = _cfit3.getFit();
+ double xc = circle.x0();
+ double yc = circle.y0();
+ double rc = Math.sqrt(xc*xc + yc*yc);
+ double rcurv = circle.radius();
+
+ // Find the point of closest approach
+ double x0 = xc * (1. - rcurv / rc);
+ double y0 = yc * (1. - rcurv / rc);
+
+ // Find the x-y arc lengths to the hits and the smallest arc length
+ double phi0 = Math.atan2(y0-yc, x0-xc);
+ double[] dphi = new double[3];
+ double dphimin = 999.;
+
+ for (int i=0; i<3; i++) {
+ // Find the angle between the hits and the DCA under the assumption that |dphi| < pi
+ dphi[i] = Math.atan2(p[i][1]-yc, p[i][0]-xc) - phi0;
+ if (dphi[i] > Math.PI) dphi[i] -= twopi;
+ if (dphi[i] < -Math.PI) dphi[i] += twopi;
+// System.out.println("dphi = "+dphi[i]+" for hit "+i);
+ if (Math.abs(dphi[i]) < Math.abs(dphimin)) dphimin = dphi[i];
+ }
+
+ // Use the hit closest to the DCA to determine the circle "direction"
+ boolean cw = dphimin < 0.;
+
+ // Find the arc lengths to the hits
+ double[] s = new double[3];
+ for (int i=0; i<3; i++) {
+
+ // Arc set to be positive if they have the same sign as dphimin
+ if (cw) s[i] = -dphi[i] * rcurv;
+ else s[i] = dphi[i] * rcurv;
+
+ // Treat the case where a point has dphi opposite in sign to dphimin as an incoming looper hit
+ if (s[i] < 0.) s[i] += twopi * rcurv;
+ }
+
+ // Order the arc lengths and z info by increasing arc length
+ for (int i=0; i<2; i++) {
+ for (int j=i+1; j<3; j++) {
+ if (s[j] < s[i]) {
+ double temp = s[i];
+ s[i] = s[j];
+ s[j] = temp;
+ temp = z[i];
+ z[i] = z[j];
+ z[j] = temp;
+ }
+ }
+ }
+
+ // Predict the middle z and see if it is consistent with the measurements
+ double slope = (z[2] - z[0]) / (s[2] - s[0]);
+ double z0 = z[0] - s[0] * slope;
+ double zpred = z0 + s[1] * slope;
+ if (Math.abs(zpred - z[1]) > dztot) return false;
+
+ // If we haven't already done the circle checks, do them now
+ if (zfirst) {
+ if (!TwoPointCircleCheck(hit1, hit3, null)) return false;
+ if (!TwoPointCircleCheck(hit2, hit3, null)) return false;
+ }
+
+ // Passed all checks - success!
+ return true;
+ }
+
+
private double dphimax(double r1, double r2) {
// Order the two radii
@@ -403,7 +492,8 @@
private void CorrectHitPosition(HelicalTrackHit hit, SeedCandidate seed) {
if (hit instanceof HelicalTrackCross) {
HelicalTrackCross cross = (HelicalTrackCross) hit;
- HelicalTrackFit helix = seed.getHelix();
+ HelicalTrackFit helix = null;
+ if (seed != null) helix = seed.getHelix();
if (helix != null) {
TrackDirection trkdir = HelixUtils.CalculateTrackDirection(helix, helix.PathMap().get(hit));
cross.setTrackDirection(trkdir, helix.covariance());
lcsim/src/org/lcsim/recon/tracking/seedtracker
diff -u -r1.6 -r1.7
--- HelixFitter.java 19 Aug 2009 22:12:27 -0000 1.6
+++ HelixFitter.java 6 Nov 2009 00:29:40 -0000 1.7
@@ -7,7 +7,6 @@
package org.lcsim.recon.tracking.seedtracker;
-import java.util.HashMap;
import java.util.List;
import java.util.Map;
@@ -18,12 +17,10 @@
import org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
import org.lcsim.fit.helicaltrack.HelixUtils;
-import org.lcsim.fit.helicaltrack.MultipleScatter;
import org.lcsim.fit.helicaltrack.TrackDirection;
import org.lcsim.fit.line.SlopeInterceptLineFit;
import org.lcsim.recon.tracking.seedtracker.diagnostic.ISeedTrackerDiagnostics;
import org.lcsim.fit.zsegment.ZSegmentFit;
-import org.omg.PortableServer._ServantActivatorStub;
/**
*
@@ -42,6 +39,7 @@
private ZSegmentFit _zsegmentfit;
private FitStatus _status;
private ISeedTrackerDiagnostics _diag = null;
+
/**
* Creates a new instance of HelixFitter
*/
lcsim/src/org/lcsim/recon/tracking/seedtracker
diff -u -r1.14 -r1.15
--- SeedTrackFinder.java 3 Sep 2009 17:52:29 -0000 1.14
+++ SeedTrackFinder.java 6 Nov 2009 00:29:40 -0000 1.15
@@ -65,53 +65,39 @@
// Instantiate the fast hit checker
FastCheck checker = new FastCheck(strategy, bfield, _diag);
- // Get the SeedLayers for this strategy
- List<SeedLayer> seedlayerlist = strategy.getLayers(SeedLayer.SeedType.Seed);
- if (seedlayerlist.size() != 3)
- throw new RuntimeException("Illegal Strategy "+strategy.getName()+": Number of Seed Layers is not 3");
-
// Find the valid sector combinations
- SeedSectoring ss = new SeedSectoring(_hitmanager, strategy, bfield, seedlayerlist);
+ SeedSectoring ss = new SeedSectoring(_hitmanager, strategy, bfield);
List<List<Sector>> sslist = ss.SeedSectors();
// Loop over the first seed layer
for (List<Sector> slist : sslist) {
for (HelicalTrackHit hit1 : slist.get(0).Hits()) {
- // Construct a seed with the current hit in the first seed layer
- SeedCandidate seed0 = new SeedCandidate(strategy, bfield);
- seed0.addHit(hit1);
-
// Loop over the second seed layer and check that we have a hit pair consistent with our strategy
for (HelicalTrackHit hit2 : slist.get(1).Hits()) {
- // Construct a new seed adding in the hit from the second seed layer
- SeedCandidate seed1 = new SeedCandidate(seed0);
- seed1.addHit(hit2);
-
- // Check if the pair of hits is consistent with current strategy
- if (!checker.CheckHitPair(hit1, hit2, seed1)) {
- if (_diag != null) _diag.fireCheckHitPairFailed(hit1, hit2, seed1);
+ // Check if the pair of hits is consistent with the current strategy
+ if (!checker.TwoPointCircleCheck(hit1, hit2, null)) {
+ if (_diag != null) _diag.fireCheckHitPairFailed(hit1, hit2);
continue;
}
// Loop over the third seed layer and check that we have a hit triplet consistent with our strategy
for (HelicalTrackHit hit3 : slist.get(2).Hits()) {
- SeedCandidate seed = new SeedCandidate(seed1);
- seed.addHit(hit3);
-
- // Check if the new hit is consistent with the other hits
- if (!checker.CheckHitPair(hit1, hit3, seed)) {
- if (_diag != null) _diag.fireCheckHitPairFailed(hit1, hit3, seed);
- continue;
- }
- if (!checker.CheckHitPair(hit2, hit3, seed)) {
- if (_diag != null) _diag.fireCheckHitPairFailed(hit2, hit3, seed);
+ // Check if the triplet of hits is consistent with the current strategy
+ if (!checker.ThreePointHelixCheck(hit1, hit2, hit3)) {
+ if (_diag != null) _diag.fireCheckHitTripletFailed(hit1, hit2, hit3);
continue;
}
- // Found a seed - if it's a true seed, add the MC Particle to those that were seeded
+ // Form a seed candidate from the seed hits
+ SeedCandidate seed = new SeedCandidate(strategy, bfield);
+ seed.addHit(hit1);
+ seed.addHit(hit2);
+ seed.addHit(hit3);
+
+ // If it's a true seed, add the MC Particle to those that were seeded
if (_diag != null) {
if (seed.isTrueSeed()) _seededmcp.addAll(seed.getMCParticles());
}
@@ -150,6 +136,7 @@
// Done with track finding for this strategy
if (_diag != null)
_diag.fireFinderDone(_trackseeds, _seededmcp);
+
return _trackseeds.size() > 0;
}
