View Javadoc

1   package org.lcsim.recon.tracking.seedtracker;
2   
3   import org.lcsim.constants.Constants;
4   import org.lcsim.fit.threepointcircle.CircleFit;
5   import org.lcsim.fit.helicaltrack.HelicalTrack2DHit;
6   import org.lcsim.fit.helicaltrack.HelicalTrack3DHit;
7   import org.lcsim.fit.helicaltrack.HelicalTrackCross;
8   import org.lcsim.fit.helicaltrack.HelicalTrackFit;
9   import org.lcsim.fit.helicaltrack.HelicalTrackHit;
10  import org.lcsim.fit.threepointcircle.ThreePointCircleFitter;
11  import org.lcsim.fit.twopointcircle.TwoPointCircleFit;
12  import org.lcsim.fit.twopointcircle.TwoPointCircleFitter;
13  import org.lcsim.fit.twopointcircle.TwoPointLineFit;
14  import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
15  import org.lcsim.recon.tracking.seedtracker.diagnostic.ISeedTrackerDiagnostics;
16  
17  /**
18   *
19   * @author Richard Partridge
20   */
21  public class FastCheck {
22  
23      private SeedStrategy _strategy;
24      private double _bfield;
25      private double _RMin;
26      private double _dMax;
27      private double _z0Max;
28      private double _nsig;
29      private TwoPointCircleFitter _cfit2;
30      private ThreePointCircleFitter _cfit3;
31      private static double twopi = 2. * Math.PI;
32      private double _eps = 1.0e-6;
33      private ISeedTrackerDiagnostics _diag;
34      private boolean _skipchecks = false;
35      private boolean _doSectorBinCheck = false;
36      private SectorManager _sectorManager = null;
37  
38      public FastCheck(SeedStrategy strategy, double bfield, ISeedTrackerDiagnostics diag) {
39          _strategy = strategy;
40          _bfield = bfield;
41          _diag = diag;
42  
43          //  Calculate the minimum radius of curvature, maximum DCA and Maximum z0
44          _RMin = strategy.getMinPT() / (Constants.fieldConversion * bfield);
45          _dMax = strategy.getMaxDCA();
46          _z0Max = strategy.getMaxZ0();
47          _nsig = Math.sqrt(strategy.getMaxChisq());
48  
49          //  Instantiate the two point circle fitter for this minimum radius, maximum DCA
50          _cfit2 = new TwoPointCircleFitter(_RMin);
51  
52          //  Instantiate the three point circle fitter
53          _cfit3 = new ThreePointCircleFitter();
54      }
55      
56      public void setDoSectorBinCheck(SectorManager sectorManager) {
57          _doSectorBinCheck = true;
58          _sectorManager = sectorManager;
59      }
60  
61      public boolean CheckHitSeed(HelicalTrackHit hit, SeedCandidate seed) {
62  
63          if (_skipchecks) return true;
64  
65          //  Check the hit against each hit in the seed
66          for (HelicalTrackHit hit2 : seed.getHits()) {
67              if (!TwoPointCircleCheck(hit, hit2, seed)) return false;
68              if (this._doSectorBinCheck) {
69                  if (!zSectorCheck(hit,hit2)) {
70                      return false;
71                  }
72              }
73          }
74  
75          return true;
76      }
77  
78      public boolean CheckSector(SeedCandidate seed, Sector sector) {
79  
80          if (_skipchecks) return true;
81  
82          //  Get limits on r, phi, and z for hits in this sector
83          double rmin = sector.rmin();
84          double rmax = sector.rmax();
85          double phimin = sector.phimin();
86          double phimax = sector.phimax();
87          double zmin = sector.zmin();
88          double zmax = sector.zmax();
89  
90          //  Calculate the midpoint and half the span in phi for this layer
91          double midphisec = (phimin + phimax) / 2.;
92          double dphisec = 0.5 * (phimax - phimin);
93  
94          //  Check each hit for compatibility with this sector
95          for (HelicalTrackHit hit : seed.getHits()) {
96  
97              //  Adjust the hit position for stereo hits
98              CorrectHitPosition(hit, seed);
99  
100             //  Calculate the max track angle change between the hit and sector layer
101             double dphitrk1 = dphimax(hit.r(), rmin);
102             double dphitrk2 = dphimax(hit.r(), rmax);
103             double dphitrk = Math.max(dphitrk1, dphitrk2);
104 
105             //  Calculate the phi dev between the hit and midpoint of the sector
106             double dphi = phidif(hit.phi(), midphisec);
107 
108             //  The maximum dphi is the sum of the track bend and half the sector span
109             double dphimx = dphitrk + dphisec;
110             if (dphi > dphimx) return false;
111 
112             double smin1 = smin(rmin);
113             double smax1 = smax(rmax);
114             double r = hit.r();
115             double smin2 = smin(r);
116             double smax2 = smax(r);
117 
118             //  Get the z limits for the hit
119             double zlen = 0.;
120             if (hit instanceof HelicalTrack2DHit) {
121                 zlen = ((HelicalTrack2DHit) hit).zlen();
122             }
123             double zmin2 = hit.z() - 0.5 * zlen;
124             double zmax2 = zmin2 + zlen;
125 
126             //  Check the z0 limits
127             boolean zOK = checkz0(smin1, smax1, zmin, zmax, smin2, smax2, zmin2, zmax2);
128             
129             if(!zOK) return false;
130             
131             boolean zSectorOK = true;
132         
133             if(_doSectorBinCheck) {                
134                 zSectorOK = zSectorCheck(hit,sector);
135             }
136             
137             if(!zSectorOK) return false;
138             
139         }
140         return true;
141     }
142 
143     public boolean CheckSectorPair(Sector s1, Sector s2) {
144 
145         if (_skipchecks) return true;
146 
147         //  Calculate the maximum change in azimuth
148         double dphi1 = dphimax(s1.rmin(), s2.rmax());
149         double dphi2 = dphimax(s1.rmax(), s2.rmin());
150 
151         //  Calculate the angular difference between the midpoints of the 2 sectors
152         double mid1 = (s1.phimax() + s1.phimin()) / 2.0;
153         double mid2 = (s2.phimax() + s2.phimin()) / 2.0;
154         double dmid = phidif(mid1, mid2);
155 
156         //  Calculate the half widths of the 2 sectors
157         double wid1 = s1.phimax() - mid1;
158         double wid2 = s2.phimax() - mid2;
159 
160         //  Check that the sectors are compatible in the bend coordinate
161         boolean phiOK;
162 
163         phiOK = dmid < dphi1 + wid1 + wid2;
164         if (!phiOK) phiOK = dmid < dphi2 + wid1 + wid2;
165         if (!phiOK) return false;
166 
167         // Get the minimum and maximum path lengths
168         double s1min = smin(s1.rmin());
169         double s2min = smin(s2.rmin());
170         double s1max = smax(s1.rmax());
171         double s2max = smax(s2.rmax());
172 
173         //  Get the minimum and maximum z's
174         double z1min = s1.zmin();
175         double z2min = s2.zmin();
176         double z1max = s1.zmax();
177         double z2max = s2.zmax();
178 
179         //  Check that the sectors are compatible in the non-bend coordinate
180         boolean zOK = checkz0(s1min, s1max, z1min, z1max, s2min, s2max, z2min, z2max);
181         
182         if (!zOK) return false;
183   
184         boolean zSectorOK = true;
185         
186         if(_doSectorBinCheck) {
187             zSectorOK = zSectorCheck(s1,s2);
188         }
189 
190         return zSectorOK;
191     }
192 
193     public boolean TwoPointCircleCheck(HelicalTrackHit hit1, HelicalTrackHit hit2, SeedCandidate seed) {
194         if (_skipchecks) return true;
195 
196         //  Initialize the hit coordinates for an unknown track direction
197         CorrectHitPosition(hit1, seed);
198         CorrectHitPosition(hit2, seed);
199 
200         //  Check that hits are outside the maximum DCA
201         if (hit1.r() < _dMax || hit2.r() < _dMax) return false;
202 
203         //  Try to find a circle passing through the 2 hits and the maximum DCA
204         boolean success = false;
205         try
206         {
207             success = _cfit2.FitCircle(hit1, hit2, _dMax);
208         }
209         catch(Exception x){}
210 
211         //  Check for success
212         if (!success) return false;
213         
214         //  Initialize the minimum/maximum arc lengths
215         double s1min = 1.0e99;
216         double s1max = -1.0e99;
217         double s2min = 1.0e99;
218         double s2max = -1.0e99;
219 
220         //  Loop over the circle fits and find the min/max arc lengths
221         for (TwoPointCircleFit fit : _cfit2.getCircleFits()) {
222             double s1 = fit.s1();
223             double s2 = fit.s2();
224             if (s1 < s1min) s1min = s1;
225             if (s1 > s1max) s1max = s1;
226             if (s2 < s2min) s2min = s2;
227             if (s2 > s2max) s2max = s2;
228         }
229 
230         //  If we are consistent with a straight-line fit, update the minimum s1, s2
231         TwoPointLineFit lfit = _cfit2.getLineFit();
232         if (lfit != null) {
233 
234             //  Find the distance from the DCA to the maximum DCA circle
235             double x0 = lfit.x0();
236             double y0 = lfit.y0();
237             double s0 = 0.;
238             double s0sq = _dMax*_dMax - (x0*x0 + y0*y0);
239             if (s0sq > _eps*_eps) s0 = Math.sqrt(s0sq);
240 
241             //  Update the minimum arc length to the distance from the DCA to the hit
242             s1min = lfit.s1() - s0;
243             s2min = lfit.s2() - s0;
244         }
245 
246         //  Calculate the allowed variation in hit r and z (not 1 sigma errors!)
247         double dr1 = Math.max(_nsig * hit1.dr(), _dMax);
248         double dr2 = Math.max(_nsig * hit2.dr(), _dMax);
249         double dz1 = dz(hit1);
250         double dz2 = dz(hit2);
251 
252         //  Now check for consistent hits in the s-z plane
253         //  First expand z ranges by hit z uncertainty
254         double z1min = hit1.z() - dz1;
255         double z1max = hit1.z() + dz1;
256         double z2min = hit2.z() - dz2;
257         double z2max = hit2.z() + dz2;
258 
259         //  Expand s ranges by hit r uncertainty (r ~ s for r << R_curvature)
260         s1min = Math.max(0., s1min - dr1);
261         s1max = s1max + dr1;
262         s2min = Math.max(0., s2min - dr2);
263         s2max = s2max + dr2;
264 
265         //  Check the z0 limits using the min/max path lengths
266         boolean zOK = checkz0(s1min, s1max, z1min, z1max, s2min, s2max, z2min, z2max);
267 
268         if(!zOK) return false;
269             
270         boolean zSectorOK = true;
271         
272         if(_doSectorBinCheck) {
273             zSectorOK = zSectorCheck(hit1,hit2);
274         }
275 
276         //  Done!
277         return zSectorOK;
278     }
279 
280     public boolean ThreePointHelixCheck(HelicalTrackHit hit1, HelicalTrackHit hit2, HelicalTrackHit hit3) {
281 
282         if (_skipchecks) return true;
283 
284         //  Setup for a 3 point circle fit
285         double p[][] = new double[3][2];
286         double[] pos;
287         double z[] = new double[3];
288         double dztot = 0.;
289         int indx;
290         HelicalTrackHit hit;
291         boolean zfirst = true;
292 
293         //  While not terribly elegant, code for speed
294         //  Use calls that give uncorrected position and error
295         //  Get the relevant variables for hit 1
296         indx = 0;
297         hit = hit1;
298         pos = hit.getPosition();
299         p[indx][0] = pos[0];
300         p[indx][1] = pos[1];
301         z[indx] = pos[2];
302         if (hit.BarrelEndcapFlag() == BarrelEndcapFlag.BARREL) {
303             if (hit instanceof HelicalTrack3DHit) dztot += _nsig * ((HelicalTrack3DHit) hit).dz();
304             else {
305                 zfirst = false;
306                 if (hit instanceof HelicalTrack2DHit) dztot += ((HelicalTrack2DHit) hit).zlen() / 2.;
307                 else dztot += _nsig * Math.sqrt(hit.getCovMatrix()[5]);
308             }
309         } else {
310             dztot += hit.dr() * Math.abs(pos[2]) / Math.sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
311         }
312 
313          //  Get the relevant variables for hit 2
314         indx = 1;
315         hit = hit2;
316         pos = hit.getPosition();
317         p[indx][0] = pos[0];
318         p[indx][1] = pos[1];
319         z[indx] = pos[2];
320         if (hit.BarrelEndcapFlag() == BarrelEndcapFlag.BARREL) {
321             if (hit instanceof HelicalTrack3DHit) dztot += _nsig * ((HelicalTrack3DHit) hit).dz();
322             else {
323                 zfirst = false;
324                 if (hit instanceof HelicalTrack2DHit) dztot += ((HelicalTrack2DHit) hit).zlen() / 2.;
325                 else dztot += _nsig * Math.sqrt(hit.getCovMatrix()[5]);
326             }
327         } else {
328             dztot += hit.dr() * Math.abs(pos[2]) / Math.sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
329         }
330 
331         //  Get the relevant variables for hit 3
332         indx = 2;
333         hit = hit3;
334         pos = hit.getPosition();
335         p[indx][0] = pos[0];
336         p[indx][1] = pos[1];
337         z[indx] = pos[2];
338         if (hit.BarrelEndcapFlag() == BarrelEndcapFlag.BARREL) {
339             if (hit instanceof HelicalTrack3DHit) dztot += _nsig * ((HelicalTrack3DHit) hit).dz();
340             else {
341                 zfirst = false;
342                 if (hit instanceof HelicalTrack2DHit) dztot += ((HelicalTrack2DHit) hit).zlen() / 2.;
343                 else dztot += _nsig * Math.sqrt(hit.getCovMatrix()[5]);
344             }
345         } else {
346             dztot += hit.dr() * Math.abs(pos[2]) / Math.sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
347         }
348 
349         //  Add multiple scattering error here - for now, just set it to 1 mm
350         dztot += 1.;
351 
352         //  Unless the three hits are all pixel hits, do the circle checks first
353         if (!zfirst) {
354             if (!TwoPointCircleCheck(hit1, hit3, null)) return false;
355             if (!TwoPointCircleCheck(hit2, hit3, null)) return false;
356         }
357 
358         //  Do the 3 point circle fit and check for success
359         boolean success = _cfit3.fit(p[0], p[1], p[2]);
360         if (!success) return false;
361 
362         //  Retrieve the circle parameters
363         CircleFit circle = _cfit3.getFit();
364         double xc = circle.x0();
365         double yc = circle.y0();
366         double rc = Math.sqrt(xc*xc + yc*yc);
367         double rcurv = circle.radius();
368 
369         //  Find the point of closest approach
370         double x0 = xc * (1. - rcurv / rc);
371         double y0 = yc * (1. - rcurv / rc);
372 
373         //  Find the x-y arc lengths to the hits and the smallest arc length
374         double phi0 = Math.atan2(y0-yc, x0-xc);
375         double[] dphi = new double[3];
376         double dphimin = 999.;
377 
378         for (int i=0; i<3; i++) {
379             //  Find the angle between the hits and the DCA under the assumption that |dphi| < pi
380             dphi[i] = Math.atan2(p[i][1]-yc, p[i][0]-xc) - phi0;
381             if (dphi[i] > Math.PI) dphi[i] -= twopi;
382             if (dphi[i] < -Math.PI) dphi[i] += twopi;
383             if (Math.abs(dphi[i]) < Math.abs(dphimin)) dphimin = dphi[i];
384             }
385 
386         //  Use the hit closest to the DCA to determine the circle "direction"
387         boolean cw = dphimin < 0.;
388 
389         //  Find the arc lengths to the hits
390         double[] s = new double[3];
391         for (int i=0; i<3; i++) {
392 
393             //  Arc set to be positive if they have the same sign as dphimin
394             if (cw) s[i] = -dphi[i] * rcurv;
395             else s[i] = dphi[i] * rcurv;
396 
397             //  Treat the case where a point has dphi opposite in sign to dphimin as an incoming looper hit
398             if (s[i] < 0.) s[i] += twopi * rcurv;
399             }
400 
401         //  Order the arc lengths and z info by increasing arc length
402         for (int i=0; i<2; i++) {
403             for (int j=i+1; j<3; j++) {
404                 if (s[j] < s[i]) {
405                     double temp = s[i];
406                     s[i] = s[j];
407                     s[j] = temp;
408                     temp = z[i];
409                     z[i] = z[j];
410                     z[j] = temp;
411                 }
412             }
413         }
414 
415         //  Predict the middle z and see if it is consistent with the measurements
416         double slope = (z[2] - z[0]) / (s[2] - s[0]);
417         double z0 = z[0] - s[0] * slope;
418         double zpred = z0 + s[1] * slope;
419         if (Math.abs(zpred - z[1]) > dztot) return false;
420 
421         //  If we haven't already done the circle checks, do them now
422         if (zfirst) {
423             if (!TwoPointCircleCheck(hit1, hit3, null)) return false;
424             if (!TwoPointCircleCheck(hit2, hit3, null)) return false;
425         }
426         
427         //  Passed all checks - success!
428         return true;
429     }
430 
431 
432     private double dphimax(double r1, double r2) {
433 
434         //  Order the two radii
435         double rmin = r1;
436         double rmax = r2;
437         if (rmin > rmax) {
438             rmin = r2;
439             rmax = r1;
440         }
441 
442         //  Don't let the maximum DCA exceed rmin
443         double d = _dMax;
444         if (d > rmin) d = rmin = _eps;
445 
446         //  First try for a circle that fits between the DCA and rmax...
447         double R = (rmax + _dMax) / 2.0 + _eps;
448         //  ...but don't let the circle have a radius smaller than the minimum allowed
449         if (R < _RMin) R = _RMin;
450 
451         //  Calculate the maximum deviations from a straight-line track
452         double cth1 = (d * d + rmin * rmin - 2 * R * d) / (2 * rmin * (R - d));
453         double cth2 = (d * d + rmax * rmax - 2 * R * d) / (2 * rmax * (R - d));
454         if (Math.abs(cth1) > 1 || Math.abs(cth2) > 1) return 0.;
455         return phidif(Math.acos(cth1), Math.acos(cth2));
456     }
457 
458     private double smin(double r) {
459         return r - _dMax;
460     }
461 
462     private double smax(double r) {
463         double R = _RMin;
464         if (r > _RMin) R = r;
465         return 2.0 * R * Math.asin((r + _dMax) / (2.0 * R));
466     }
467 
468     private boolean checkz0(double s1min, double s1max, double zmin1, double zmax1,
469             double s2min, double s2max, double zmin2, double zmax2) {
470 
471         double z0[] = new double[2];
472         double z1[] = new double[4];
473         double z2[] = new double[2];
474         double s1[] = new double[4];
475         double s2[] = new double[2];
476 
477         //  Set limits on z0
478         z0[0] = -_z0Max;
479         z0[1] = _z0Max;
480 
481         //  Set corners of allowed region for s1, z1
482         z1[0] = zmin1;
483         z1[1] = zmin1;
484         z1[2] = zmax1;
485         z1[3] = zmax1;
486         s1[0] = s1min;
487         s1[1] = s1max;
488         s1[2] = s1min;
489         s1[3] = s1max;
490 
491         //  Set limits on z2, s2
492         z2[0] = zmin2;
493         z2[1] = zmax2;
494         s2[0] = s2min;
495         s2[1] = s2max;
496 
497         //  Initialize min/max of s, z at point 2
498         double zmax = -1.0e10;
499         double zmin = 1.0e10;
500         double smax = -1.0e10;
501         double smin = 1.0e10;
502         
503         //  Loop over z0 limits
504         for (int i=0; i<2; i++) {
505 
506             //  Loop over corners of s1, z1
507             for (int j=0; j<4; j++) {
508 
509                 //  Calculate slope of line in s-z space from z0 limit to point 1 corner
510                 double slope = (z1[j] - z0[i]) / s1[j];
511 
512                 //  Loop over limits on z2, s2
513                 for (int k=0; k<2; k++) {
514 
515                     //  Calculate extrapolation of s-z line to the point 2 limit
516                     double z = z0[i] + s2[k] * slope;
517                     double s = (z2[k] - z0[i]) / slope;
518 
519                     //  Find the min/max values of the extrapolated s, z at point 2
520                     if (z > zmax) zmax = z;
521                     if (z < zmin) zmin = z;
522                     if (s > smax) smax = s;
523                     if (s < smin) smin = s;
524                 }
525             }
526         }
527 
528         //  Check to see if the extrapolated points are consistent with measurements
529         boolean checkz0 = (zmin2 <= zmax && zmax2 >= zmin) || (s2min <= smax && s2max >= smin);
530 
531         return checkz0;
532     }
533 
534     private double phidif(double phi1, double phi2) {
535 
536         //  Find the magnitude of the phi difference
537         double phidif = Math.abs(phi2 - phi1);
538 
539         //  Properly wrap around two pi to always give the smaller phi dif
540         if (phidif > Math.PI) phidif = 2. * Math.PI - phidif;
541         return phidif;
542     }
543 
544     private void CorrectHitPosition(HelicalTrackHit hit, SeedCandidate seed) {
545         if (hit instanceof HelicalTrackCross) {
546             HelicalTrackCross cross = (HelicalTrackCross) hit;
547             HelicalTrackFit helix = null;
548             if (seed != null) helix = seed.getHelix();
549             cross.setTrackDirection(helix);
550         }
551     }
552 
553     private double dz(HelicalTrackHit hit) {
554 
555         //  Axial strip hits: use half strip length
556         if (hit instanceof HelicalTrack2DHit) {
557             return 0.5 * ((HelicalTrack2DHit) hit).zlen();
558 
559         //  Otherwise use the z error
560         } else {
561             return _nsig * Math.sqrt(hit.getCorrectedCovMatrix().diagonal(2));
562         }
563     }
564 
565     private boolean zSectorCheck(Sector s1, Sector s2) {
566         return s1.zSector()==s2.zSector();
567     }
568 
569     private boolean zSectorCheck(HelicalTrackHit hit, Sector sector) {
570         int zSector = sector.zSector();
571         int zBin = this._sectorManager.ZBin(hit);
572         return zBin==zSector;
573     }
574 
575     private boolean zSectorCheck(HelicalTrackHit hit, HelicalTrackHit hit2) {
576         int zBin = this._sectorManager.ZBin(hit);
577         int zBin2 = this._sectorManager.ZBin(hit2);
578         return zBin==zBin2;
579     }
580 
581 }