lcsim/src/org/lcsim/recon/tracking/seedtracker
diff -u -r1.8 -r1.9
--- ConfirmerExtender.java 7 Feb 2009 03:56:06 -0000 1.8
+++ ConfirmerExtender.java 19 Feb 2009 22:00:34 -0000 1.9
@@ -100,7 +100,6 @@
// Loop over the confirmation/extension layers
for (SeedLayer lyr : lyrs) {
-// System.out.println("New Layer: "+lyr.LayerID());
// Create a list for seeds that are confirmed/extended using the current layer
List<SeedCandidate> newlist = new ArrayList<SeedCandidate>();
@@ -119,31 +118,19 @@
double oldcirclechisq = seedin.getHelix().chisq()[0];
double chisqbest = 1.e6;
- // Use the first seed hit for the fast check of additional hits
-// HelicalTrackHit checkhit = seedin.getHits().get(0);
-
-// System.out.println("Trial seed with "+seedin.getHits().size()+" hits");
-// System.out.println("Sectors with hits: "+sectorlist.size());
-
// Create a list of hits to check for this layer / seed
List<HelicalTrackHit> hitlist = new ArrayList<HelicalTrackHit>();
// Loop over the sectors with hits
for (Sector sector : sectorlist) {
-// System.out.println("About to check a sector");
// Check that this candidate is compatible with the sector
-// if (!checker.CheckSector(seedin, sector)) continue;
-
-// System.out.println("Sector checked out OK");
-// System.out.println("Found a sector to check: "+
-// sector.Identifier()+" with "+sector.Hits().size()+" hits");
+ if (!checker.CheckSector(seedin, sector)) continue;
// Add the hits for this confirmation/extension sector to the hitlist
hitlist.addAll(sector.Hits());
}
-// System.out.println("Hits to check: "+hitlist.size());
SortHits comp = new SortHits(seedin.getHelix());
Collections.sort(hitlist, comp);
@@ -154,9 +141,8 @@
// For each hit in this confirmation/extension layer, make a test seed including the new hit
for (HelicalTrackHit hit : hitlist) {
-// System.out.println("Checking a hit for viability");
// Check that this hit is potentially viable
-// if (!checker.CheckHitSeed(hit, seedin)) continue;
+ if (!checker.CheckHitSeed(hit, seedin)) continue;
SeedCandidate test = new SeedCandidate(seedin);
test.addHit(hit);
@@ -164,7 +150,6 @@
// See if we have a successful fit
boolean success = fitter.FitCandidate(test, strategy);
-// System.out.println("Fit completed with status: "+success);
if (success) {
// Attach the fit to the test seed and add it to the list
@@ -180,12 +165,10 @@
if (diag != null) {
diag.fireConfirmerExtenderFitNoSuccess(task, seedin, hit, fitter, true);
}
-// System.out.println("Fit Status: "+fitter.getFitStatus());
+
+ // Stop checking hits in this layer if circle chisq increase is too big
if (fitter.getFitStatus() != FitStatus.CircleFitFailed) {
double circlechisq = fitter.getCircleFit().chisq();
-
-// System.out.println("Delta chisq: "+(circlechisq-oldcirclechisq));
-
if (circlechisq > oldcirclechisq + strategy.getMaxChisq())
break;
}
@@ -195,7 +178,6 @@
// Finished checking hits to add to this seed candidate
// If all fit tries for this layer are potentially bad hits, include the starting seed in the list
if (chisqbest - oldchisq > strategy.getBadHitChisq()) {
-// System.out.println("Keeping original fit with "+seedin.getHits().size()+" hits");
newlist.add(seedin);
// if (diag != null) {
// diag.fireConfirmerExtenderAllHitsPotentiallyBad(task, seedin, hitlist, chisqbest, oldchisq);
@@ -208,7 +190,6 @@
// Finished looping over the seeds we are trying to confirm/extend with this layer
// Use the new list of seeds as the working listinput to the next layer search
seedlist = newlist;
-// System.out.println("New list with "+seedlist.size()+" seeds");
if (diag != null) {
diag.fireConfirmerExtenderLayerDone(task, n, seedlist);
}
@@ -221,7 +202,6 @@
// Check for seeds with sufficient numbers of confirmed/extended hits
for (SeedCandidate candidate : seedlist) {
int hits = candidate.getHits().size();
-// System.out.println("Final candidate with "+hits+" hits");
if (hits >= minhits) result.add(candidate);
}
return result.size() > 0;
lcsim/src/org/lcsim/recon/tracking/seedtracker
diff -u -r1.5 -r1.6
--- FastCheck.java 7 Feb 2009 03:56:06 -0000 1.5
+++ FastCheck.java 19 Feb 2009 22:00:34 -0000 1.6
@@ -1,8 +1,11 @@
package org.lcsim.recon.tracking.seedtracker;
+import java.util.List;
import org.lcsim.constants.Constants;
+import org.lcsim.event.MCParticle;
import org.lcsim.fit.helicaltrack.HelicalTrack2DHit;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
/**
*
@@ -14,13 +17,7 @@
private double _RMin;
private double _dMax;
private double _z0Max;
- private double _r1;
- private double _phi1;
- private double _dphimx1;
- private double _zmin1;
- private double _zmax1;
- private double _s1min;
- private double _s1max;
+ private double _eps = 1.0e-6;
public FastCheck(SeedStrategy strategy, double bfield) {
// Calculate the minimum radius of curvature, maximum DCA and Maximum z0
@@ -31,54 +28,67 @@
public boolean CheckHitPair(HelicalTrackHit hit1, HelicalTrackHit hit2) {
- // If we have a cache miss for hit1, calculate the cached hit 1 variables
- if (hit1 != _hit1cache) {
- _r1 = hit1.r();
- _phi1 = hit1.phi();
- _dphimx1 = dphimx(_r1);
- double zlen1 = 0.;
- if (hit1 instanceof HelicalTrack2DHit) {
- zlen1 = ((HelicalTrack2DHit) hit1).zlen();
- }
- _zmin1 = hit1.z() - 0.5 * zlen1;
- _zmax1 = _zmin1 + zlen1;
- _s1min = smin(_r1);
- _s1max = smax(_r1);
- _hit1cache = hit1;
- }
-
- // Get the polar coordinates for hit 2
+ // 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 from a radial track given the cuts on radius of curvature and DCA
- double dphimx2 = dphimx(r2);
+ // 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);
- double dphimx = phidif(dphimx2, _dphimx1);
+ double dphi = phidif(phi2, phi1);
// Check that the difference in azimuthal angle is less than the difference in maximum phi deviations
- if (dphi > dphimx) return false;
+ boolean phiOK = dphi <= dphimx;
+
+ if (!phiOK) {
+// if (CheckMC(hit1, hit2)) {
+// 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 hit2
+ // 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(_s1min, _s1max, _zmin1, _zmax1, smin(r2), smax(r2), zmin2, zmax2);
+ boolean zOK = checkz0(smin(r1), smax(r1), zmin1, zmax1, smin(r2), smax(r2), zmin2, zmax2);
+// if (!zOK) {
+// if (CheckMC(hit1, hit2)) {
+// 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);
+// }
+// }
return zOK;
}
public boolean CheckHitSeed(HelicalTrackHit hit, SeedCandidate seed) {
// Check the hit against each hit in the seed
- for (HelicalTrackHit hit2: seed.getHits()) {
+ for (HelicalTrackHit hit2 : seed.getHits()) {
if (!CheckHitPair(hit, hit2)) return false;
}
return true;
@@ -94,10 +104,6 @@
double zmin = sector.zmin();
double zmax = sector.zmax();
- // Calculate the maximum phi deviation from a straight-line track for this layer
- double dphimxsec1 = dphimx(rmax);
- double dphimxsec2 = dphimx(rmin);
-
// Calculate the midpoint and half the span in phi for this layer
double midphisec = (phimin + phimax) / 2.;
double dphisec = 0.5 * (phimax - phimin);
@@ -106,11 +112,10 @@
for (HelicalTrackHit hit : seed.getHits()) {
// Calculate the max track angle change between the hit and sector layer
- double dphimxhit = dphimx(hit.r());
- double dphitrk1 = phidif(dphimxsec1, dphimxhit);
- double dphitrk2 = phidif(dphimxsec2, dphimxhit);
+ double dphitrk1 = dphimax(hit.r(), rmin);
+ double dphitrk2 = dphimax(hit.r(), rmax);
double dphitrk = Math.max(dphitrk1, dphitrk2);
-
+
// Calculate the phi dev between the hit and midpoint of the sector
double dphi = phidif(hit.phi(), midphisec);
@@ -134,7 +139,6 @@
// Check the z0 limits
boolean zOK = checkz0(smin1, smax1, zmin, zmax, smin2, smax2, zmin2, zmax2);
-// System.out.println("zOK: "+zOK);
if (!zOK) return false;
}
return true;
@@ -143,12 +147,8 @@
public boolean CheckSectorPair(Sector s1, Sector s2) {
// Calculate the maximum change in azimuth
- double dphi1rmin = dphimx(s1.rmin());
- double dphi2rmin = dphimx(s2.rmin());
- double dphi1rmax = dphimx(s1.rmax());
- double dphi2rmax = dphimx(s2.rmax());
- double dphi1 = phidif(dphi1rmax, dphi2rmin);
- double dphi2 = phidif(dphi1rmin, dphi2rmax);
+ double dphi1 = dphimax(s1.rmin(), s2.rmax());
+ double dphi2 = dphimax(s1.rmax(), s2.rmin());
// Calculate the angular difference between the midpoints of the 2 sectors
double mid1 = (s1.phimax() + s1.phimin()) / 2.0;
@@ -161,11 +161,10 @@
// Check that the sectors are compatible in the bend coordinate
boolean phiOK;
-// System.out.println("Checking phi");
+
phiOK = dmid < dphi1 + wid1 + wid2;
if (!phiOK) phiOK = dmid < dphi2 + wid1 + wid2;
if (!phiOK) return false;
-// System.out.println("Phi OK");
// Get the minimum and maximum path lengths
double s1min = smin(s1.rmin());
@@ -182,15 +181,33 @@
// Check that the sectors are compatible in the non-bend coordinate
boolean zOK = checkz0(s1min, s1max, z1min, z1max, s2min, s2max, z2min, z2max);
-// System.out.println("sector pair check OK: "+zOK);
-
return zOK;
}
- private double dphimx(double r) {
- double R = _RMin;
- if (r > _RMin) R = r;
- return Math.asin((2 * R * _dMax - _dMax * _dMax - r * r) / (2 * r * (R - _dMax)));
+ private double dphimax(double r1, double r2) {
+
+ // Order the two radii
+ double rmin = r1;
+ double rmax = r2;
+ if (rmin > rmax) {
+ rmin = r2;
+ rmax = r1;
+ }
+
+ // Don't let the maximum DCA exceed rmin
+ double d = _dMax;
+ if (d > rmin) d = rmin = _eps;
+
+ // First try for a circle that fits between the DCA and rmax...
+ double R = (rmax + _dMax) / 2.0 + _eps;
+ // ...but don't let the circle have a radius smaller than the minimum allowed
+ if (R < _RMin) R = _RMin;
+
+ // Calculate the maximum deviations from a straight-line track
+ double cth1 = (d * d + rmin * rmin - 2 * R * d) / (2 * rmin * (R - d));
+ double cth2 = (d * d + rmax * rmax - 2 * R * d) / (2 * rmax * (R - d));
+ if (Math.abs(cth1) > 1 || Math.abs(cth2) > 1) return 0.;
+ return phidif(Math.acos(cth1), Math.acos(cth2));
}
private double smin(double r) {
@@ -198,25 +215,28 @@
}
private double smax(double r) {
- return 2.0 * _RMin * Math.asin((r + _dMax) / (2.0 * _RMin));
+ double R = _RMin;
+ if (r > _RMin) R = r;
+ return 2.0 * R * Math.asin((r + _dMax) / (2.0 * R));
}
private boolean checkz0(double s1min, double s1max, double zmin1, double zmax1,
double s2min, double s2max, double zmin2, double zmax2) {
// Find the limits on z0 for smin and smax
- double zlim[] = new double[4];
- zlim[0] = (zmin1 * s2min - zmax2 * s1min) / (s2min - s1min);
- zlim[1] = (zmax1 * s2min - zmin2 * s1min) / (s2min - s1min);
- zlim[2] = (zmin1 * s2max - zmax2 * s1max) / (s2max - s1max);
- zlim[3] = (zmax1 * s2max - zmin2 * s1max) / (s2max - s1max);
-
- // Find the largest and smallest z0 limits
- double zmax = -99999.;
- double zmin = 99999.;
- for (int i = 0; i<4; i++) {
- if (zlim[i]>zmax) zmax = zlim[i];
- if (zlim[i]<zmin) zmin = zlim[i];
+ double zlim[] = new double[6];
+ zlim[0] = (zmin1 * s2min - zmax2 * s1max) / (s2min - s1max);
+ zlim[1] = (zmax1 * s2max - zmin2 * s1min) / (s2max - s1min);
+ zlim[2] = (zmin1 * s2max - zmax2 * s1min) / (s2max - s1min);
+ zlim[3] = (zmax1 * s2min - zmin2 * s1max) / (s2min - s1max);
+ zlim[4] = (zmin1 * s2max - zmin2 * s1max) / (s2max - s1max);
+ zlim[5] = (zmax1 * s2min - zmax2 * s1min) / (s2min - s1min);
+
+ double zmin = zlim[0];
+ double zmax = zlim[0];
+ for (int i=1; i<6; i++) {
+ if (zlim[i] < zmin) zmin = zlim[i];
+ if (zlim[i] > zmax) zmax = zlim[i];
}
// System.out.println("s1min: "+s1min+" s2: "+s2min+" zmin1: "+zmin1+
@@ -239,4 +259,23 @@
return phidif;
}
+ private boolean CheckMC(HelicalTrackHit hit1, HelicalTrackHit hit2) {
+ List<MCParticle> mclist1 = hit1.getMCParticles();
+ List<MCParticle> mclist2 = hit2.getMCParticles();
+ for (MCParticle mc1 : mclist1) {
+ HelixParamCalculator helix = new HelixParamCalculator(mc1, 5.0);
+ if (helix.getMCTransverseMomentum() < 0.2) continue;
+ if (Math.abs(helix.getDCA()) > 10.) continue;
+ if (Math.abs(helix.getZ0()) > 10.) continue;
+ if (mclist2.contains(mc1)) {
+ System.out.println("pt: "+helix.getMCTransverseMomentum()+" p: "+helix.getMCMomentum());
+ System.out.println(" MC z0 origin: "+mc1.getOriginZ()+" p: "+mc1.getMomentum().magnitude());
+ double s1 = (hit1.getPosition()[2] - helix.getZ0()) / helix.getSlopeSZPlane();
+ double s2 = (hit2.getPosition()[2] - helix.getZ0()) / helix.getSlopeSZPlane();
+ System.out.println("s1: "+s1+" s2: "+s2);
+ return true;
+ }
+ }
+ return false;
+ }
}
\ No newline at end of file