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
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
44 _RMin = strategy.getMinPT() / (Constants.fieldConversion * bfield);
45 _dMax = strategy.getMaxDCA();
46 _z0Max = strategy.getMaxZ0();
47 _nsig = Math.sqrt(strategy.getMaxChisq());
48
49
50 _cfit2 = new TwoPointCircleFitter(_RMin);
51
52
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
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
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
91 double midphisec = (phimin + phimax) / 2.;
92 double dphisec = 0.5 * (phimax - phimin);
93
94
95 for (HelicalTrackHit hit : seed.getHits()) {
96
97
98 CorrectHitPosition(hit, seed);
99
100
101 double dphitrk1 = dphimax(hit.r(), rmin);
102 double dphitrk2 = dphimax(hit.r(), rmax);
103 double dphitrk = Math.max(dphitrk1, dphitrk2);
104
105
106 double dphi = phidif(hit.phi(), midphisec);
107
108
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
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
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
148 double dphi1 = dphimax(s1.rmin(), s2.rmax());
149 double dphi2 = dphimax(s1.rmax(), s2.rmin());
150
151
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
157 double wid1 = s1.phimax() - mid1;
158 double wid2 = s2.phimax() - mid2;
159
160
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
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
174 double z1min = s1.zmin();
175 double z2min = s2.zmin();
176 double z1max = s1.zmax();
177 double z2max = s2.zmax();
178
179
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
197 CorrectHitPosition(hit1, seed);
198 CorrectHitPosition(hit2, seed);
199
200
201 if (hit1.r() < _dMax || hit2.r() < _dMax) return false;
202
203
204 boolean success = false;
205 try
206 {
207 success = _cfit2.FitCircle(hit1, hit2, _dMax);
208 }
209 catch(Exception x){}
210
211
212 if (!success) return false;
213
214
215 double s1min = 1.0e99;
216 double s1max = -1.0e99;
217 double s2min = 1.0e99;
218 double s2max = -1.0e99;
219
220
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
231 TwoPointLineFit lfit = _cfit2.getLineFit();
232 if (lfit != null) {
233
234
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
242 s1min = lfit.s1() - s0;
243 s2min = lfit.s2() - s0;
244 }
245
246
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
253
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
260 s1min = Math.max(0., s1min - dr1);
261 s1max = s1max + dr1;
262 s2min = Math.max(0., s2min - dr2);
263 s2max = s2max + dr2;
264
265
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
277 return zSectorOK;
278 }
279
280 public boolean ThreePointHelixCheck(HelicalTrackHit hit1, HelicalTrackHit hit2, HelicalTrackHit hit3) {
281
282 if (_skipchecks) return true;
283
284
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
294
295
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
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
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
350 dztot += 1.;
351
352
353 if (!zfirst) {
354 if (!TwoPointCircleCheck(hit1, hit3, null)) return false;
355 if (!TwoPointCircleCheck(hit2, hit3, null)) return false;
356 }
357
358
359 boolean success = _cfit3.fit(p[0], p[1], p[2]);
360 if (!success) return false;
361
362
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
370 double x0 = xc * (1. - rcurv / rc);
371 double y0 = yc * (1. - rcurv / rc);
372
373
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
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
387 boolean cw = dphimin < 0.;
388
389
390 double[] s = new double[3];
391 for (int i=0; i<3; i++) {
392
393
394 if (cw) s[i] = -dphi[i] * rcurv;
395 else s[i] = dphi[i] * rcurv;
396
397
398 if (s[i] < 0.) s[i] += twopi * rcurv;
399 }
400
401
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
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
422 if (zfirst) {
423 if (!TwoPointCircleCheck(hit1, hit3, null)) return false;
424 if (!TwoPointCircleCheck(hit2, hit3, null)) return false;
425 }
426
427
428 return true;
429 }
430
431
432 private double dphimax(double r1, double r2) {
433
434
435 double rmin = r1;
436 double rmax = r2;
437 if (rmin > rmax) {
438 rmin = r2;
439 rmax = r1;
440 }
441
442
443 double d = _dMax;
444 if (d > rmin) d = rmin = _eps;
445
446
447 double R = (rmax + _dMax) / 2.0 + _eps;
448
449 if (R < _RMin) R = _RMin;
450
451
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
478 z0[0] = -_z0Max;
479 z0[1] = _z0Max;
480
481
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
492 z2[0] = zmin2;
493 z2[1] = zmax2;
494 s2[0] = s2min;
495 s2[1] = s2max;
496
497
498 double zmax = -1.0e10;
499 double zmin = 1.0e10;
500 double smax = -1.0e10;
501 double smin = 1.0e10;
502
503
504 for (int i=0; i<2; i++) {
505
506
507 for (int j=0; j<4; j++) {
508
509
510 double slope = (z1[j] - z0[i]) / s1[j];
511
512
513 for (int k=0; k<2; k++) {
514
515
516 double z = z0[i] + s2[k] * slope;
517 double s = (z2[k] - z0[i]) / slope;
518
519
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
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
537 double phidif = Math.abs(phi2 - phi1);
538
539
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
556 if (hit instanceof HelicalTrack2DHit) {
557 return 0.5 * ((HelicalTrack2DHit) hit).zlen();
558
559
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 }