00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include <algorithm>
00016 #include <string.h>
00017 #include "CLHEP/String/Strings.h"
00018
00019 #include "TrkReco/TMDC.h"
00020 #include "TrkReco/TMDCUtil.h"
00021 #include "TrkReco/TTrackManager.h"
00022 #include "TrkReco/TConformalFinder.h"
00023 #include "TrkReco/TConformalFinder0.h"
00024 #include "TrkReco/THistogram.h"
00025 #include "TrkReco/TMDCTsf.h"
00026 #include "TrkReco/TCircle.h"
00027 #include "TrkReco/TTrack.h"
00028 #include "TrkReco/TMLine.h"
00029 #include "TrkReco/TTrackHEP.h"
00030 #include "TrkReco/TMDCWireHitMC.h"
00031
00032 #ifdef TRKRECO_DEBUG
00033 unsigned TConformalFinder::_stage = ConformalOutside;
00034 #endif
00035 #ifdef TRKRECO_WINDOW
00036 #include "TrkReco/TWindow.h"
00037 #endif
00038
00039
00040
00041 static const float PI2 = 2. * M_PI;
00042 static const float WIDTH[] = {PI2 / 40,
00043 PI2 / 64,
00044 PI2 / 76,
00045 PI2 / 100,
00046 PI2 / 128,
00047 PI2 / 160,
00048 PI2 / 176,
00049 PI2 / 208,
00050 PI2 / 240,
00051 PI2 / 256,
00052 PI2 / 288};
00053
00054 std::string
00055 TConformalFinder::version(void) const {
00056 return "3.03";
00057 }
00058
00059 TConformalFinder::TConformalFinder(unsigned fastFinder,
00060 unsigned slowFinder,
00061 unsigned perfectSegmentFinding,
00062 float maxSigma,
00063 float maxSigmaStereo,
00064 float salvageLevel,
00065 unsigned minNLinksForSegment,
00066 unsigned minNCoreLinks,
00067 unsigned minNSegments,
00068 unsigned salvageLoadWidth,
00069 unsigned stereoMode,
00070 unsigned stereoLoadWidth,
00071 float szSegmentDistance,
00072 float szLinkDistance,
00073 unsigned fittingFlag)
00074 : _fastFinder(fastFinder),
00075 _slowFinder(slowFinder),
00076 _doT0Reset(false),
00077 _perfectSegmentFinding(perfectSegmentFinding),
00078 _segmentSeparation(4),
00079 _minNLinksForSegment(minNLinksForSegment),
00080 _minNLinksForSegmentInRefine(2),
00081 _maxNLinksForSegment(8),
00082 _maxWidthForSegment(4),
00083 _minNCoreLinks(minNCoreLinks),
00084 _minNSegments(minNSegments),
00085 _minUsedFractionSlow2D(0.5),
00086 _salvageLoadWidth(salvageLoadWidth),
00087 _stereoMode(stereoMode),
00088 _stereoLoadWidth(stereoLoadWidth),
00089 _appendLoad(2),
00090 _maxSigma2(maxSigma),
00091 _T0ResetDone(false),
00092 _builder("conformal builder",
00093 maxSigma,
00094 maxSigmaStereo,
00095 salvageLevel,
00096 szSegmentDistance,
00097 szLinkDistance,
00098 fittingFlag),
00099 _s(0)
00100 #ifdef TRKRECO_WINDOW
00101 , _rphiWindow("rphi window")
00102 #endif
00103 {
00104 _linkMaxDistance[0] = 0.02;
00105 _linkMaxDistance[1] = 0.025;
00106 _linkMaxDistance[2] = 0.025;
00107 _linkMinDirAngle[0] = 0.97;
00108 _linkMinDirAngle[1] = 0.95;
00109 _linkMinDirAngle[2] = 0.95;
00110
00111
00112
00113 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
00114 if(scmgn!=StatusCode::SUCCESS) {
00115 std::cout<< __FILE__<<"Unable to open Magnetic field service"<<std::endl;
00116 }
00117 }
00118
00119 TConformalFinder::~TConformalFinder() {
00120 }
00121
00122 void
00123 TConformalFinder::dump(const std::string & msg, const std::string & pre) const {
00124 if (msg.find("state") != std::string::npos) {
00125 std::cout << pre;
00126 TFinderBase::dump(msg);
00127 std::cout << pre;
00128 std::cout << "#axial=" << _hits[0].length();
00129 std::cout << ",#stereo=" << _hits[1].length();
00130 }
00131 if (msg.find("summary") != std::string::npos || msg.find("detail") != std::string::npos) {
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 }
00170 }
00171
00172 void
00173 TConformalFinder::clear(void) {
00174 for (unsigned i = 0; i < 3; i++) {
00175 if (i == 2)
00176 HepAListDeleteAll(_allHits[i]);
00177 else
00178 _allHits[i].removeAll();
00179 _hits[i].removeAll();
00180 _unused[i].removeAll();
00181 }
00182 for (unsigned i = 0; i < 2; i++) {
00183 for (unsigned j = 0; j < 6; j++) {
00184 #ifdef TRKRECO_DEBUG_DETAIL
00185 cout<<"_allSegments length = "<<_allSegments[i][j].length()<<" _allUnused length = "<<_allUnused[i][j].length()<<endl;
00186 #endif
00187 HepAListDeleteAll(_allSegments[i][j]);
00188 _allUnused[i][j].removeAll();
00189 }
00190 }
00191 _2DTracks.removeAll();
00192 _3DTracks.removeAll();
00193 }
00194
00195 void
00196 TConformalFinder::selectGoodHits(void) {
00197 for (unsigned i = 0; i < 2; i++) {
00198 unsigned n = _allHits[i].length();
00199 for (unsigned j = 0; j < n; j++) {
00200 unsigned state = _allHits[i][j]->hit()->state();
00201
00202
00203
00204 _hits[i].append(_allHits[i][j]);
00205
00206 _unused[i].append(_allHits[i][j]);
00207 }
00208 }
00209 _hits[2].append(_hits[0]);
00210 _hits[2].append(_hits[1]);
00211 _unused[2].append(_unused[0]);
00212 _unused[2].append(_unused[1]);
00213
00214 #ifdef TRKRECO_WINDOW
00215 _rphiWindow.clear();
00216 _rphiWindow.skip(false);
00217 _rphiWindow.skipAllWindow(false);
00218 _rphiWindow.append(_allHits[2]);
00219 _rphiWindow.append(_hits[2], leda_pink);
00220
00221
00222 #endif
00223 }
00224
00225 void
00226 TConformalFinder::findSegments(void) {
00227
00228
00229 AList<TMLink> links[11];
00230 unsigned n = _hits[2].length();
00231 for (unsigned i = 0; i < n; i++) {
00232 TMLink & l = * _hits[2][i];
00233 links[l.wire()->superLayerId()].append(l);
00234 }
00235
00236
00237 THistogram * hist[11];
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 hist[0] = new THistogram(56);
00251 hist[1] = new THistogram(80);
00252 hist[2] = new THistogram(88);
00253 hist[3] = new THistogram(112);
00254 hist[4] = new THistogram(140);
00255 hist[5] = new THistogram(160);
00256 hist[6] = new THistogram(176);
00257 hist[7] = new THistogram(208);
00258 hist[8] = new THistogram(240);
00259 hist[9] = new THistogram(256);
00260 hist[10] = new THistogram(288);
00261
00262 for (unsigned i = 0; i < 11; i++) {
00263
00264
00265 unsigned superLayerId;
00266 unsigned id;
00267 switch(i){
00268 case 0: superLayerId = 0; id = 1; break;
00269 case 1: superLayerId = 1; id = 1; break;
00270 case 2: superLayerId = 0; id = 0; break;
00271 case 3: superLayerId = 1; id = 0; break;
00272 case 4: superLayerId = 2; id = 0; break;
00273 case 5: superLayerId = 2; id = 1; break;
00274 case 6: superLayerId = 3; id = 1; break;
00275 case 7: superLayerId = 4; id = 1; break;
00276 case 8: superLayerId = 5; id = 1; break;
00277 case 9: superLayerId = 3; id = 0; break;
00278 case 10: superLayerId = 4; id = 0; break;
00279 default: break;
00280 }
00281
00282 hist[i]->fillPhi(links[i]);
00283
00284
00285 AList<TSegment> tmp = hist[i]->segments();
00286
00287
00288 unsigned n = tmp.length();
00289 if (id == 0) {
00290 AList<TSegment> bad;
00291 for (unsigned j = 0; j < n; j++) {
00292 TSegment * s = tmp[j];
00293 if (s->links().length() < _minNLinksForSegment)
00294 bad.append(s);
00295 else if (s->links().length() > _maxNLinksForSegment)
00296 bad.append(s);
00297
00298
00299 }
00300 tmp.remove(bad);
00301 for (unsigned j = 0; j < bad.length(); j++) {
00302 _unused[id].append(bad[j]->links());
00303 _unused[2].append(bad[j]->links());
00304 }
00305 HepAListDeleteAll(bad);
00306 n = tmp.length();
00307 }
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 _allSegments[id][superLayerId].append(tmp);
00329 _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
00330 delete hist[i];
00331 }
00332
00333 #ifdef TRKRECO_DEBUG_DETAIL
00334 std::cout << "... segment finding finished" << std::endl;
00335 #endif
00336 }
00337
00338 void
00339 TConformalFinder::linkSegments(unsigned level) {
00340
00341 unsigned superLayer = 5;
00342 while (--superLayer) {
00343 AList<TSegment> & segments = _allUnused[0][superLayer];
00344 unsigned n = segments.length();
00345
00346 for (unsigned i = 0; i < n; i++) {
00347 TSegment & base = * segments[i];
00348 #ifdef TRKRECO_DEBUG
00349 if (base.tracks().length()) {
00350 std::cout << "TConformalFinder::linkSegments !!! segment has ";
00351 std::cout << "an assignment to track(s)" << std::endl;
00352 }
00353 #endif
00354 base.innerLinks().removeAll();
00355 base.outerLinks().removeAll();
00356 }
00357 }
00358
00359
00360 superLayer = 5;
00361 while (--superLayer) {
00362 AList<TSegment> & segments = _allUnused[0][superLayer];
00363 AList<TSegment> & targets = _allUnused[0][superLayer - 1];
00364 AList<TSegment> & targets2 = _allUnused[0][superLayer - 2];
00365 unsigned n = segments.length();
00366 for (unsigned i = 0; i < n; i++) {
00367 TSegment & base = * segments[i];
00368 #ifdef TRKRECO_DEBUG
00369 if (base.tracks().length()) {
00370 std::cout << "TConformalFinder::linkSegments !!! segment has ";
00371 std::cout << "an assignment to track(s)" << std::endl;
00372 }
00373 #endif
00374
00375 const HepPoint3D & p = base.position();
00376 HepVector3D v = base.direction();
00377
00378 if (base.outerLinks().length() == 1)
00379 v = p - OuterMostUniqueLink(base)->position();
00380
00381
00382
00383
00384
00385 AList<TSegment> & candidates = base.innerLinks();
00386 TSegment * best = link(base, p, v, targets, candidates, level);
00387 if ((best == NULL) && (level > 0) && (superLayer > 1)) {
00388 best = link(base, p, v, targets2, candidates, level);
00389 }
00390 if (best) candidates.insert(best);
00391
00392
00393
00394 unsigned m = candidates.length();
00395 for (unsigned j = 0; j < m; j++)
00396 candidates[j]->outerLinks().append(base);
00397 }
00398 }
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436 #ifdef TRKRECO_DEBUG
00437 superLayer = 5;
00438 while(--superLayer){
00439 cout<<"SuperLayer: "<<superLayer<<endl;
00440 AList<TSegment> & segments = _allUnused[0][superLayer];
00441 unsigned n = segments.length();
00442 for (unsigned i = 0; i < n; i++) {
00443 TSegment & base = * segments[i];
00444 cout<<"Layer: "<<base.links()[0]->hit()->wire()->layerId()
00445 <<" Local: "<<base.links()[0]->hit()->wire()->localId()
00446 <<" innerSeg: "<<base.innerLinks().length()<<endl;
00447 }
00448 }
00449 #endif
00450
00451 #ifdef TRKRECO_WINDOW
00452
00453 displayStatus("Conf::linkSegments ... results");
00454 _rphiWindow.wait();
00455 #endif
00456 }
00457
00458 void
00459 TConformalFinder::resolveSegments(AList<TTrack> & trackCandidates) const {
00460 #ifdef TRKRECO_DEBUG
00461 std::cout << "... resolving assignments" << std::endl;
00462 #endif
00463
00464
00465 AList<TSegment> multi;
00466 unsigned n = trackCandidates.length();
00467 for (unsigned i = 0; i < n; i++) {
00468 TTrack & t = * trackCandidates[i];
00469 AList<TSegment> & segments = t.segments();
00470 unsigned nS = segments.length();
00471 for (unsigned j = 0; j < nS; j++) {
00472 if (segments[j]->tracks().length() > 1)
00473 multi.append(segments[j]);
00474 }
00475 }
00476 multi.purge();
00477
00478
00479 n = multi.length();
00480 for (unsigned i = 0; i < n; i++) {
00481 TSegment & s = * multi[i];
00482 const AList<TTrack> & tracks = s.tracks();
00483 unsigned nT = tracks.length();
00484
00485
00486 AList<TMLink> multiLinks;
00487 const AList<TMLink> & links = s.links();
00488 unsigned nL = links.length();
00489 for (unsigned j = 0; j < nL; j++) {
00490 TMLink & l = * links[i];
00491 bool overlap = false;
00492 for (unsigned k = 0; k < nT; k++) {
00493 TTrack & t = * tracks[k];
00494 if (t.links().hasMember(l))
00495 overlap = true;
00496 }
00497 multiLinks.append(l);
00498 }
00499 multiLinks.purge();
00500
00501 #ifdef TRKRECO_DEBUG
00502 std::cout << " segment : ";
00503 s.dump("hits sort flag");
00504 std::cout << " # of assigned tracks = " << nT << std::endl;
00505 std::cout << " # of overlap TMLinks = " << multiLinks.length();
00506 std::cout << std::endl;
00507 #endif
00508
00509 nL = multiLinks.length();
00510 for (unsigned j = 0; j < nL; j++) {
00511 TMLink & l = * multiLinks[j];
00512
00513 float bestDiff = 999999999.;
00514 TTrack * best = NULL;
00515 for (unsigned k = 0; k < nT; k++) {
00516 TTrack & t = * tracks[k];
00517 t.approach(l);
00518 float distance = (l.positionOnWire() - l.positionOnTrack())
00519 .mag();
00520 float diff = fabs(distance - l.hit()->drift());
00521 if (diff < bestDiff) {
00522 bestDiff = diff;
00523 best = & t;
00524 }
00525 }
00526
00527 for (unsigned k = 0; k < nT; k++) {
00528 TTrack & t = * tracks[k];
00529 if (& t == best) continue;
00530 t.remove(l);
00531 }
00532
00533 #ifdef TRKRECO_DEBUG
00534 {
00535 std::cout << " " << l.wire()->name() << " -> ";
00536 std::cout << best->name() << std::endl;
00537 }
00538 #endif
00539 }
00540 }
00541 }
00542
00543 AList<TSegment>
00544 TConformalFinder::removeBadSegments(TTrack & t) const {
00545 #ifdef TRKRECO_DEBUG
00546 std::cout << "... removing bad segments : # used seg = ";
00547 std::cout << t.segments().length() << std::endl;
00548 #endif
00549 #ifdef TRKRECO_DEBUG_DETAIL
00550 for (unsigned i = 0; i < t.segments().length(); i++)
00551 t.segments()[i]->dump("hits sort flag", " ");
00552 #endif
00553
00554 const AList<TSegment> & segments = t.segments();
00555 AList<TSegment> bads;
00556 unsigned used = 0;
00557 TSegment * innerMost;
00558 unsigned n = segments.length();
00559 for (unsigned i = 0; i < n; i++) {
00560 TSegment & s = * segments[i];
00561 AList<TMLink> links = Links(s, t);
00562 unsigned nLinks = links.length();
00563
00564 unsigned nCores = Cores(links).length();
00565 unsigned nCoresSegment = s.nCores();
00566
00567
00568 if ((nCores < _minNCoreLinks) && (nCores < ((nCoresSegment + 1) / 2))){
00569 bads.append(s);
00570 continue;
00571 }
00572
00573 unsigned id = segments[i]->superLayerId();
00574 used |= (1 << id);
00575 if (! id) innerMost = & s;
00576 }
00577
00578
00579
00580
00581
00582
00583
00584 if (! bads.length()) return bads;
00585
00586 #ifdef TRKRECO_DEBUG
00587 std::cout << " bad segments : # = " << bads.length() << std::endl;
00588 #endif
00589
00590 n = bads.length();
00591 for (unsigned i = 0; i < n; i++) {
00592 TSegment & s = * bads[i];
00593
00594 #ifdef TRKRECO_DEBUG
00595 AList<TMLink> links = Links(s, t);
00596 unsigned nLinks = links.length();
00597 unsigned nCores = Cores(links).length();
00598 std::cout << " # used links = " << nLinks;
00599 std::cout << ", # used cores = " << nCores << std::endl;
00600 s.dump("hits sort flag", " ");
00601 #endif
00602
00603 s.tracks().remove(t);
00604 t.segments().remove(s);
00605 t.remove(s.links());
00606 }
00607
00608
00609 t.fit();
00610
00611 #ifdef TRKRECO_DEBUG
00612 std::cout << "... refitting" << std::endl;
00613 t.dump("detail", "2nd> ");
00614 #endif
00615
00616 return bads;
00617 }
00618
00619 TSegment *
00620 TConformalFinder::link(const TSegment & base,
00621 const HepPoint3D & p,
00622 const HepVector3D & v,
00623 const AList<TSegment> & candidates,
00624 AList<TSegment> & alternatives,
00625 unsigned level) const {
00626
00627 #ifdef TRKRECO_DEBUG
00628 std::cout << " link : base = " << std::endl;
00629 base.dump("vector hits mc sort", "-> ");
00630 std::cout << " p=" << p << ", v=" << v.unit() << std::endl;
00631 #endif
00632
00633
00634
00635
00636
00637
00638 unsigned n = candidates.length();
00639 float maxAngle = -999.;
00640 float minDistance = _linkMaxDistance[level];
00641 float maxDirAngle = _linkMinDirAngle[level];
00642
00643 TSegment * best = NULL;
00644 for (unsigned j = 0; j < n; j++) {
00645 TSegment & current = * candidates[j];
00646
00647 #ifdef TRKRECO_DEBUG_DETAIL
00648 current.dump("vector hits mc sort", " ");
00649 std::cout << " dist,dirAngle,angle=" << current.distance(p, v);
00650 std::cout << "," << ((current.position() - p).unit()).dot(v.unit());
00651 std::cout << "," << v.dot(current.direction()) << std::endl;
00652 #endif
00653 float distance = current.distance(p, v);
00654
00655
00656
00657 if (distance > _linkMaxDistance[level]) continue;
00658
00659 HepVector3D dir = (current.position() - p);
00660 if (dir.x() > M_PI) dir.setX(dir.x() - PI2);
00661 else if (dir.x() < - M_PI) dir.setX(PI2 + dir.x());
00662
00663 float dirAngle = dir.unit().dot(v.unit());
00664
00665 if (dirAngle < _linkMinDirAngle[level]) continue;
00666
00667
00668
00669
00670
00671 if (dirAngle > maxDirAngle) {
00672 maxDirAngle = dirAngle;
00673 best = & current;
00674 }
00675 alternatives.append(current);
00676 }
00677
00678 #ifdef TRKRECO_DEBUG
00679 std::cout << " # of best + alternatives = " << alternatives.length() << std::endl;
00680 if (best) {
00681 std::cout << " Best is ";
00682 best->dump("hits");
00683 }
00684 #endif
00685
00686 alternatives.remove(best);
00687 if (best) return best;
00688 return NULL;
00689 }
00690
00691 void
00692 TConformalFinder::deleteTrack(TTrack & t) const {
00693 const AList<TSegment> & segments = t.segments();
00694 unsigned n = segments.length();
00695 for (unsigned i = 0; i < n; i++) {
00696 TSegment & s = * segments[i];
00697 s.tracks().remove(t);
00698 }
00699
00700 delete & t;
00701 }
00702
00703 void
00704 TConformalFinder::removeUsedSegments(const AList<TTrack> & tracks) {
00705 unsigned n = tracks.length();
00706 for (unsigned i = 0; i < n; i++) {
00707 TTrack & t = * tracks[i];
00708 AList<TSegment> & segments = t.segments();
00709 AList<TSegment> toBeRemoved;
00710 unsigned nS = segments.length();
00711 for (unsigned j = 0; j < nS; j++) {
00712 TSegment & s = * segments[j];
00713 unsigned sId = s.superLayerId();
00714
00715
00716 unsigned as;
00717 unsigned id;
00718 switch (sId) {
00719 case 0: as = 1; id = 0;break;
00720 case 1: as = 1; id = 1;break;
00721 case 2: as = 0; id = 0;break;
00722 case 3: as = 0; id = 1;break;
00723 case 4: as = 0; id = 2;break;
00724 case 5: as = 1; id = 2;break;
00725 case 6: as = 1; id = 3;break;
00726 case 7: as = 1; id = 4;break;
00727 case 8: as = 1; id = 5;break;
00728 case 9: as = 0; id = 3;break;
00729 case 10: as = 0; id = 4;break;
00730 default: break;
00731 }
00732
00733
00734 AList<TMLink> links = Links(s, t);
00735 if (links.length() == 0) {
00736 s.tracks().remove(t);
00737 toBeRemoved.append(s);
00738 #ifdef TRKRECO_DEBUG
00739 std::cout << "!!! why this happends???" << std::endl;
00740 std::cout << " no used link" << std::endl;
00741 #endif
00742 continue;
00743 }
00744
00745
00746 _allUnused[as][id].remove(s);
00747
00748
00749 const AList<TSegment> & outers = s.outerLinks();
00750 unsigned nO = outers.length();
00751
00752
00753 for (unsigned i = 0; i < nO; i++) {
00754 outers[i]->innerLinks().remove(s);
00755
00756 }
00757
00758
00759 AList<TMLink> unused = s.links();
00760 unused.remove(links);
00761 s.remove(unused);
00762 unsigned nUnused = unused.length();
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787 }
00788 segments.remove(toBeRemoved);
00789 }
00790 }
00791
00792 void
00793 TConformalFinder::updateTLinks(AList<TTrack> & tracks) {
00794 unsigned n = tracks.length();
00795 for (unsigned i = 0; i < n; i++) {
00796 const AList<TMLink> & links = tracks[i]->links();
00797 unsigned nL = links.length();
00798 for (unsigned j = 0; j < nL; j++) {
00799 TMLink & l = * links[j];
00800 tracks[i]->approach(l);
00801 }
00802 }
00803 }
00804
00805 int
00806 TConformalFinder::doit(const AList<TMDCWireHit> & axial,
00807 const AList<TMDCWireHit> & stereo,
00808 AList<TTrack> & tracks,
00809 AList<TTrack> & tracks2D) {
00810 #ifdef TRKRECO_DEBUG
00811 _stage = ConformalInitialization;
00812 #endif
00813
00814 static bool first = true;
00815 if (first) {
00816 first = false;
00817 int size;
00818
00819
00820
00821
00822 }
00823
00824
00825 if (debugLevel()) {
00826 std::cout << name() << " ... processing"
00827 << " axial=" << axial.length()
00828 << ",stereo=" << stereo.length()
00829 << ",tracks=" << tracks.length()
00830 << std::endl;
00831 }
00832
00833
00834
00835
00836
00837 TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, axial, _allHits[0]);
00838 TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, stereo, _allHits[1]);
00839 _allHits[2].append(_allHits[0]);
00840 _allHits[2].append(_allHits[1]);
00841
00842
00843 selectGoodHits();
00844 #ifdef TRKRECO_DEBUG
00845 cout<<"axial Hits:"<<_allHits[0].length()<<" good axial hits:"<<_hits[0].length()
00846 <<endl<<"stereo Hits:"<<_allHits[1].length()<<" good setero hits:"<<_hits[1].length()
00847 <<endl;
00848 #endif
00849
00850
00851 if (_perfectSegmentFinding)
00852 findSegmentsPerfect();
00853 else
00854
00855 findSegmentsTsf();
00856
00857 #ifdef TRKRECO_DEBUG
00858 cout<<"axial Seg: "<<_allSegments[0][0].length()<<", "<<_allSegments[0][1].length()
00859 <<", "<<_allSegments[0][2].length()<<", "<<_allSegments[0][3].length()
00860 <<", "<<_allSegments[0][4].length()<<", "<<_allSegments[0][5].length()
00861 <<" stereo Seg: "<<_allSegments[1][0].length()<<", "<<_allSegments[1][1].length()
00862 <<", "<<_allSegments[1][2].length()<<", "<<_allSegments[1][3].length()
00863 <<", "<<_allSegments[1][4].length()<<", "<<_allSegments[1][5].length()<<endl;
00864
00865 for(int i=0; i<2; ++i) {
00866 if (i == 0) cout<<"Axial......"<<endl;
00867 else cout<<"Stereo......"<<endl;
00868 for(int j=0; j<6; ++j) {
00869 for(int k=0; k<_allSegments[i][j].length(); ++k){
00870 cout<<"SL(a/s): "<<j<<" seeds in Seg"<<k<<": "<<_allSegments[i][j][k]->links().length()<<endl;
00871 for(int kk=0; kk<_allSegments[i][j][k]->links().length(); ++kk)
00872 cout<<" layerId: "<<_allSegments[i][j][k]->links()[kk]->hit()->wire()->layerId()
00873 <<" localId: "<<_allSegments[i][j][k]->links()[kk]->hit()->wire()->localId()<<endl;
00874 }
00875 }
00876 }
00877 #endif
00878
00879
00880 unsigned level = 0;
00881 _T0ResetDone = false;
00882 if (_fastFinder) {
00883 linkSegments(level);
00884 fastFinding2D(level);
00885 updateTLinks(_2DTracks);
00886
00887
00888 if (_doT0Reset) {
00889 std::cout << "TConformalFinder ... T0 reset is done" << std::endl;
00890 _T0ResetDone = true;
00891 return 0;
00892 }
00893
00894 #ifdef TRKRECO_WINDOW
00895 _rphiWindow.skip(false);
00896 #endif
00897 fastFinding3D(level);
00898 updateTLinks(_2DTracks);
00899 updateTLinks(_3DTracks);
00900
00901 #ifdef TRKRECO_WINDOW
00902 _rphiWindow.skip(false);
00903 displayTracks(_2DTracks, _allUnused, "all 2D after fast level 0");
00904 displayTracks(_3DTracks, _allUnused, "all 3D after fast level 0");
00905 #endif // TRKRECO_WINDOW
00906
00907
00908
00909 level = 1;
00910 linkSegments(level);
00911 fastFinding2D(level);
00912 updateTLinks(_2DTracks);
00913
00914 #ifdef TRKRECO_WINDOW
00915 _rphiWindow.skip(false);
00916 #endif
00917 fastFinding3D(level);
00918 updateTLinks(_2DTracks);
00919 updateTLinks(_3DTracks);
00920
00921 #ifdef TRKRECO_WINDOW
00922 _rphiWindow.skip(false);
00923 displayTracks(_2DTracks, _allUnused, "all 2D after fast level 1");
00924 displayTracks(_3DTracks, _allUnused, "all 3D after fast level 1");
00925 #endif //TRKRECO_WINDOW
00926 }
00927
00928
00929 if (_slowFinder) {
00930 level = 2;
00931 linkSegments(level);
00932 slowFinding2D(level);
00933 updateTLinks(_2DTracks);
00934 #ifdef TRKRECO_WINDOW
00935 _rphiWindow.skip(false);
00936
00937 #endif
00938 fastFinding3D(level);
00939 updateTLinks(_2DTracks);
00940 updateTLinks(_3DTracks);
00941
00942 #ifdef TRKRECO_WINDOW
00943 _rphiWindow.skip(false);
00944
00945 displayTracks(_2DTracks, _allUnused, "all 2D after slow level 2");
00946 displayTracks(_3DTracks, _allUnused, "all 3D after slow level 2");
00947 #endif //TRKRECO_WINDOW
00948 }
00949
00950 #ifdef TRKRECO_DEBUG
00951 int zsltrk = _2DTracks.length();
00952 std::cout << name() << " ... # 3D tracks = " << _3DTracks.length()
00953 << ", # 2D tracks = " << _2DTracks.length() << std::endl;
00954 for (unsigned i = 0; i < zsltrk; i++){
00955
00956 cout<<"pt:"<<_2DTracks[i]->pt()<<" impact:"<<_2DTracks[i]->impact()<<endl;
00957
00958 }
00959 #endif
00960
00961
00962 TTrackManager::maskBadHits(_3DTracks, _maxSigma2);
00963
00964 tracks2D = _2DTracks;
00965 tracks = _3DTracks;
00966 for(int i=0;i<_3DTracks.length();i++){
00967 _3DTracks[i]->finder(TrackOldConformalFinder);
00968 _3DTracks[i]->setFinderType(2);
00969 }
00970 #ifdef TRKRECO_DEBUG
00971
00972
00973 for(int nn = 0; nn < tracks2D.length(); ++nn){
00974 cout<<"2D: "<<nn<<" radius: "<<tracks2D[nn]->radius()<<endl;
00975 for (int mm = 0; mm < tracks2D[nn]->links().length(); ++mm)
00976 cout<<"layer: "<<tracks2D[nn]->links()[mm]->wire()->layerId()
00977 <<" local: "<<tracks2D[nn]->links()[mm]->wire()->localId()<<endl;
00978 }
00979
00980 cout<<"unused axial Seg: "<<_allUnused[0][0].length()<<", "<<_allUnused[0][1].length()
00981 <<", "<<_allUnused[0][2].length()<<", "<<_allUnused[0][3].length()
00982 <<", "<<_allUnused[0][4].length()<<", "<<_allUnused[0][5].length()
00983 <<" unused stereo Seg: "<<_allUnused[1][0].length()<<", "<<_allUnused[1][1].length()
00984 <<", "<<_allUnused[1][2].length()<<", "<<_allUnused[1][3].length()
00985 <<", "<<_allUnused[1][4].length()<<", "<<_allUnused[1][5].length()
00986 <<endl;
00987
00988 for(int i=0; i<2; ++i) {
00989 if (i == 0) cout<<"unused axial hits: "<<endl;
00990 else cout<<"unused stereo hits: "<<endl;
00991 for(int k = 0; k < _allHits[i].length(); ++k) {
00992 if (_allHits[i][k]->hit()->state() & WireHitUsed) continue;
00993 else cout<<" layerId: "<<_allHits[i][k]->hit()->wire()->layerId()
00994 <<" localId: "<<_allHits[i][k]->hit()->wire()->localId()<<endl;
00995 }
00996 if (i == 0) cout<<"unused axial link in segs: "<<endl;
00997 else cout<<"unused stereo link in segs: "<<endl;
00998 for(int j=0; j<6; ++j) {
00999 for(int k=0; k<_allUnused[i][j].length(); ++k){
01000 cout<<"sl: "<<i<<" "<<j<<" length of seg "<<k<<": "<<_allUnused[i][j][k]->links().length()<<endl;
01001 for(int kk=0; kk<_allUnused[i][j][k]->links().length(); ++kk)
01002 cout<<" layerId: "<<_allUnused[i][j][k]->links()[kk]->hit()->wire()->layerId()
01003 <<" localId: "<<_allUnused[i][j][k]->links()[kk]->hit()->wire()->localId()<<endl;
01004 }
01005 }
01006 }
01007
01008 #endif
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036 return 0;
01037 }
01038
01039 void
01040 TConformalFinder::fastFinding3D(unsigned level) {
01041
01042 #ifdef TRKRECO_DEBUG_DETAIL
01043 std::cout << "... fast finding (3D) : level = " << level << std::endl;
01044 #endif
01045 #ifdef TRKRECO_DEBUG
01046 if (_stage == ConformalSlow2D) {
01047 _stage = ConformalSlow3D;
01048 }
01049 else {
01050 _stage = ConformalFast3DLevel0;
01051 if (level == 1)
01052 _stage = ConformalFast3DLevel1;
01053 }
01054 #endif
01055
01056 AList<TTrack> tracks3D;
01057 AList<TTrack> touched;
01058 AList<TSegment> bads;
01059 unsigned n = _2DTracks.length();
01060 for (unsigned i = 0; i < n; i++) {
01061
01062 const TTrack & t = * _2DTracks[i];
01063
01064 #ifdef TRKRECO_DEBUG
01065 cout<<"links in 2Dtracks: "<<t.links().length()<<endl;
01066 for (int ii = 0; ii < t.links().length(); ++ii) {
01067 cout<<"layerId: "<<t.links()[ii]->wire()->layerId()
01068 <<" localId: "<<t.links()[ii]->wire()->localId()<<endl;
01069 }
01070 cout<<"center: "<<t.center()<<" radius = "<<t.radius()<<endl;
01071 #endif
01072
01073 #ifdef TRKRECO_DEBUG_DETAIL
01074 std::cout << "==> fast 3D building : " << t.name() << std::endl;
01075 t.dump("hits sort flag", " 2D track hits=");
01076 #endif
01077 AList<TSegment> segments = stereoSegments(t);
01078
01079 unsigned NS = segments.length();
01080 if ( NS >= 30 ) continue;
01081
01082
01083
01084
01085
01086
01087
01088
01089 AList<TSegment> badSegments;
01090 if (_stereoMode == 2)
01091 badSegments = stereoSegmentsFromBadHits(t);
01092
01093
01094
01095
01096
01097
01098
01099 #ifdef TRKRECO_WINDOW
01100 displayStatus("Conf::fast3D ... seed");
01101 _rphiWindow.append(segments, leda_blue);
01102 _rphiWindow.append(badSegments, leda_red);
01103 _rphiWindow.oneShot(t, leda_green);
01104 #endif
01105
01106
01107
01108
01109 TTrack * s = NULL;
01110 if (_stereoMode)
01111 s = _builder.buildStereoNew(t, segments, badSegments);
01112 else
01113 s = _builder.buildStereo(t, segments);
01114 HepAListDeleteAll(badSegments);
01115
01116 if (! s) {
01117 #ifdef TRKRECO_DEBUG_DETAIL
01118 std::cout << "... 3D failure" << std::endl;
01119 #endif
01120 #ifdef TRKRECO_WINDOW
01121 displayStatus("Conf::fastd3D ... 3D failed");
01122 _rphiWindow.append(segments, leda_blue);
01123 _rphiWindow.oneShot(t, leda_red);
01124 #endif
01125 continue;
01126 }
01127
01128 if (! TTrackManager::goodTrack(* s)) {
01129 #ifdef TRKRECO_DEBUG_DETAIL
01130 std::cout << "... 3D failure (bad quality)" << std::endl;
01131 #endif
01132 #ifdef TRKRECO_WINDOW
01133 displayStatus("Conf::fastd3D ... 3D failed (bad quality)");
01134 _rphiWindow.append(segments, leda_blue);
01135 _rphiWindow.oneShot(* s, leda_red);
01136 #endif
01137
01138
01139
01140
01141
01142
01143
01144
01145 deleteTrack(* s);
01146 continue;
01147 }
01148
01149
01150 s->name(t.name() + "-3D");
01151
01152 #ifdef TRKRECO_WINDOW
01153 displayStatus("Conf::fastd3D ... 3D ok");
01154 _rphiWindow.append(segments, leda_blue);
01155 _rphiWindow.oneShot(* s, leda_green);
01156 #endif
01157
01158
01159 salvage(* s, 3, bads);
01160
01161 tracks3D.append(s);
01162 touched.append(_2DTracks[i]);
01163
01164 s->assign(WireHitConformalFinder);
01165 s->quality(0);
01166
01167
01168 static AList<TTrack> tmp;
01169 tmp.removeAll();
01170 tmp.append(s);
01171 removeUsedSegments(tmp);
01172
01173 if(tracks3D.length()>20)break;
01174 #ifdef TRKRECO_DEBUG_DETAIL
01175 std::cout << "... 3D finished : " << s->name() << std::endl;
01176 s->dump("detail", " ");
01177 #endif
01178 #ifdef TRKRECO_WINDOW
01179 displayStatus("Conf::fastd3D ... finished");
01180 _rphiWindow.oneShot(* s, leda_green);
01181 #endif
01182 }
01183
01184
01185 _3DTracks.append(tracks3D);
01186 _2DTracks.remove(tracks3D);
01187 _2DTracks.remove(touched);
01188 for (unsigned i = 0; i < touched.length(); i++) {
01189 deleteTrack(* touched[i]);
01190
01191 }
01192 }
01193
01194 #ifdef TRKRECO_WINDOW
01195 void
01196 TConformalFinder::displayStatus(const std::string & m) const {
01197 _rphiWindow.clear();
01198 _rphiWindow.text(m);
01199 _rphiWindow.append(_allHits[2], leda_pink);
01200 _rphiWindow.append(_hits[2], leda_cyan);
01201 _rphiWindow.append(_unused[2]);
01202 displayAppendSegments(_allSegments, leda_grey1);
01203 displayAppendSegments(_allUnused, leda_orange);
01204 }
01205
01206 void
01207 TConformalFinder::displayAppendSegments(const AList<TSegment> a[2][6],
01208 leda_color c) const {
01209 for (unsigned i = 0; i < 2; i++) {
01210 for (unsigned j = 0; j < 6; j++) {
01211 const AList<TSegment> & segments = a[i][j];
01212 unsigned n = segments.length();
01213 for (unsigned k = 0; k < n; k++) {
01214 _rphiWindow.append(* segments[k], c);
01215 }
01216 }
01217 }
01218 }
01219
01220 void
01221 TConformalFinder::displayTracks(const AList<TTrack> & l,
01222 const AList<TSegment> segments[2][6],
01223 const std::string & c) const {
01224 displayStatus("Conf::display ... " + c);
01225 for (unsigned i = 0; i < l.length(); i++)
01226 _rphiWindow.append(* l[i], leda_green);
01227 _rphiWindow.wait();
01228 }
01229 #endif
01230
01231 bool
01232 TConformalFinder::quality2D(TTrack & t, unsigned level) const {
01233 #ifdef TRKRECO_WINDOW
01234 std::string s = "Conf::quality2D ... level " + itostring(level);
01235 #endif
01236
01237
01238 if(fabs(t.radius()) < 15.) return false;
01239
01240
01241
01242 unsigned nSuperLayers = NSuperLayers(t.links());
01243 if (nSuperLayers < _minNSegments) {
01244 #ifdef TRKRECO_DEBUG
01245 std::cout << "... quality check : level=" << level << " : bad" << std::endl;
01246 std::cout << " reason : # segments(superlayer) =" << nSuperLayers;
01247 std::cout << " < " << _minNSegments << std::endl;
01248 #endif
01249 #ifdef TRKRECO_WINDOW
01250 s += " rejected because of few segments(superlayers)";
01251 displayStatus(s);
01252 _rphiWindow.oneShot(t, leda_red);
01253 #endif
01254 return false;
01255 }
01256 if (level == 0) {
01257 #ifdef TRKRECO_WINDOW
01258 s += " ok : # of used segments(superlayers) = ";
01259 s += itostring(nSuperLayers);
01260 s += " > " + itostring(_minNSegments);
01261 displayStatus(s);
01262 _rphiWindow.oneShot(t, leda_green);
01263 #endif
01264
01265 return true;
01266 }
01267
01268
01269 unsigned n3 = NSuperLayers(t.links(), 3);
01270
01271
01272 if (n3 != nSuperLayers) {
01273 unsigned sl = SuperLayer(t.links(), 2);
01274
01275 if (sl == 0) {
01276 #ifdef TRKRECO_DEBUG
01277 std::cout << "... quality check : level = " << level << " : bad";
01278 std::cout << std::endl;
01279 std::cout << " reason : super layer pattern(n2) = 0 " << std::endl;
01280 #endif
01281 #ifdef TRKRECO_WINDOW
01282 s += " rejected because of bad super-layer pattern(n2=0)";
01283 displayStatus(s);
01284 _rphiWindow.oneShot(t, leda_red);
01285 #endif
01286 return false;
01287 }
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321 }
01322
01323 #ifdef TRKRECO_DEBUG
01324 std::cout << "... quality check : level = " << level;
01325 std::cout << " : ok, super layer pattern = " << std::endl;
01326 std::string ptn;
01327 unsigned sl = SuperLayer(t.links(), 2);
01328 for (unsigned j = 0; j < 6; j++) {
01329 unsigned k = (sl >> (j * 2)) % 2;
01330 std::cout << k;
01331 if (k) ptn += "1";
01332 else ptn += "0";
01333 }
01334 std::cout << std::endl;
01335 #endif
01336 #ifdef TRKRECO_WINDOW
01337 s += " ok : super layer ptn = " + ptn;
01338 displayStatus(s);
01339 _rphiWindow.oneShot(t, leda_green);
01340 #endif
01341
01342 return true;
01343 }
01344
01345 void
01346 TConformalFinder::fastFinding2D(unsigned level) {
01347
01348 #ifdef TRKRECO_DEBUG_DETAIL
01349 std::cout << "... fast finding (2D) : level = " << level << std::endl;
01350 #endif
01351 #ifdef TRKRECO_DEBUG
01352 _stage = ConformalFast2DLevel0;
01353 if (level == 1)
01354 _stage = ConformalFast2DLevel1;
01355 #endif
01356
01357 AList<TSegment> * segments = & _allUnused[0][0];
01358
01359 unsigned idBase = _2DTracks.length() + _3DTracks.length();
01360 unsigned nTrial = 0;
01361 unsigned seedLayer = 5;
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372 while ((seedLayer--) > 1) {
01373
01374
01375 AList<TTrack> trackCandidates;
01376 unsigned n = segments[seedLayer].length();
01377
01378 for (unsigned i = 0; i < n; i++) {
01379 TSegment & seed = * segments[seedLayer][i];
01380 std::string name = "f" + itostring(nTrial + idBase);
01381
01382
01383
01384 #ifdef TRKRECO_DEBUG
01385 std::cout << ">>> fast 2D : super layer " << seedLayer << ", ";
01386 std::cout << i << " : ";
01387 seed.dump("link");
01388 #endif
01389
01390
01391 AList<TSegment> seeds;
01392
01393
01394 if (NUniqueLinks(seed) > 0) {
01395 seeds = UniqueLinks(seed);
01396 }
01397 else if (NMajorLinks(seed) > 0) {
01398 seeds = MajorLinks(seed);
01399 }
01400 else {
01401 continue;
01402 }
01403 seeds.append(seed);
01404
01405 #ifdef TRKRECO_DEBUG
01406 cout<<"seeds in fastFinding2D: "<<seeds.length()<<endl;
01407 for (int kk = 0; kk < seeds.length(); ++kk)
01408 cout<<"LayerId: "<<seeds[kk]->links()[0]->wire()->layerId()
01409 <<" LocalId: "<<seeds[kk]->links()[0]->wire()->localId()<<endl;
01410 #endif
01411
01412
01413 refine:
01414
01415
01416 for (int k = 0; k < seeds.length(); ++k)
01417 seeds[k]->solveLR();
01418
01419 #ifdef TRKRECO_WINDOW
01420 displayStatus("Conf::fast2D ... seed segments");
01421 _rphiWindow.oneShot(seeds, leda_green);
01422 #endif
01423
01424
01425 #ifdef TRKRECO_DEBUG
01426 std::cout << "==> fast building 2D : seed layer = " << seedLayer;
01427 std::cout << ", " << name << std::endl;
01428 #endif
01429
01430 TTrack * t = _builder.buildRphi(seeds);
01431 ++nTrial;
01432
01433
01434
01435 if (! t) continue;
01436 t->name(name);
01437 bool ok = quality2D(* t, 0);
01438
01439 if (! ok) {
01440 deleteTrack(* t);
01441 continue;
01442 }
01443
01444
01445
01446 AList<TSegment> bads = removeBadSegments(* t);
01447
01448
01449
01450 if (bads.length()) {
01451 ok = quality2D(* t, 1);
01452 if (! ok) {
01453 seeds = refineSegments(* t);
01454 if (seeds.length() >= _minNSegments) {
01455 deleteTrack(* t);
01456 goto refine;
01457 }
01458 deleteTrack(* t);
01459 continue;
01460 }
01461 }
01462
01463 salvage(* t, 1, bads);
01464
01465
01466 refineLinks(* t, 2);
01467
01468
01469
01470
01471 ok = quality2D(* t, 2);
01472 if (! ok) {
01473 deleteTrack(* t);
01474 continue;
01475 }
01476
01477
01478
01479 if (NSuperLayers(t->links()) < 5) {
01480 t = expand(* t);
01481 }
01482
01483
01484
01485 #ifdef TRKRECO_DEBUG
01486 std::cout << "... 2D finished : " << t->name() << std::endl;
01487 t->dump("hits sort flag pull", " ");
01488 #endif
01489 #ifdef TRKRECO_WINDOW
01490 displayStatus("Conf::fast2D ... finished");
01491 _rphiWindow.oneShot(* t, leda_green);
01492 #endif
01493 t->finder(TrackFastFinder);
01494 t->quality(TrackQuality2D);
01495 t->fitting(TrackFitGlobal);
01496 trackCandidates.append(t);
01497 }
01498
01499
01500 resolveSegments(trackCandidates);
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512 removeUsedSegments(trackCandidates);
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522 _2DTracks.append(trackCandidates);
01523 }
01524
01525 resolveHits(_2DTracks);
01526 }
01527
01528 const TMDCWire *
01529 TConformalFinder::conformal2Wire(const HepPoint3D & p) {
01530 static const TMDC & cdc = * TMDC::getTMDC();
01531
01532 std::cout << "p = " << p << std::endl;
01533 float r = sqrt(4. / p.mag2());
01534 float phi = p.phi();
01535 return cdc.wire(r, phi);
01536 }
01537
01538 AList<TSegment>
01539 TConformalFinder::pickUpSegmentsInConformal(float phi[12],
01540 float loadWidth,
01541 unsigned axialStereoSwitch) const {
01542 AList<TSegment> outList;
01543 HepPoint3D center(1., 1., 0.);
01544 center.setPhi(phi[0]);
01545
01546
01547 for (unsigned sl = 0; sl < 2; sl++) {
01548 if (sl == 0) {
01549 if (! (axialStereoSwitch & 1)) continue;
01550 }
01551 else {
01552 if (! (axialStereoSwitch & 2)) continue;
01553 }
01554 for (unsigned i = 0; i < 6; i++) {
01555 const AList<TSegment> & list = _allUnused[sl][i];
01556 unsigned n = list.length();
01557 for (unsigned j = 0; j < n; j++) {
01558 TSegment & s = * list[j];
01559
01560 const HepPoint3D & p = s.position();
01561 if (center.dot(p) < 0.) continue;
01562
01563 bool found = false;
01564 const AList<TMLink> & inners = s.inners();
01565 unsigned m = inners.length();
01566 for (unsigned k = 0; k < m; k++) {
01567 float dPhi = phi[i * 2 + sl] -
01568 inners[k]->position().phi();
01569 dPhi = fabs(fmod(dPhi, PI2));
01570 if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
01571 if (dPhi < WIDTH[i * 2 + sl] * loadWidth) {
01572 outList.append(s);
01573 found = true;
01574 break;
01575 }
01576 }
01577
01578 if (found) continue;
01579 const AList<TMLink> & outers = s.outers();
01580 m = outers.length();
01581 for (unsigned k = 0; k < m; k++) {
01582 float dPhi = phi[i * 2 + sl + 1] -
01583 outers[k]->position().phi();
01584 dPhi = fabs(fmod(dPhi, PI2));
01585 if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
01586 if (dPhi < WIDTH[i * 2 + sl] * loadWidth) {
01587 outList.append(s);
01588 break;
01589 }
01590 }
01591 }
01592 }
01593 }
01594
01595 return outList;
01596 }
01597
01598 AList<TMLink>
01599 TConformalFinder::pickUpLinksInConformal(float phi[12],
01600 float loadWidth,
01601 unsigned axialStereoSwitch) const {
01602 HepPoint3D center(1., 1., 0.);
01603 center.setPhi(phi[0]);
01604
01605 AList<TMLink> outList;
01606 unsigned nBad = _unused[2].length();
01607 for (unsigned i = 0; i < nBad; i++) {
01608 unsigned sl = _unused[2][i]->wire()->superLayerId();
01609 unsigned as = sl % 2;
01610 if (as == 0) {
01611 if (! (axialStereoSwitch & 1)) continue;
01612 }
01613 else {
01614 if (! (axialStereoSwitch & 2)) continue;
01615 }
01616
01617 const HepPoint3D & p = _unused[2][i]->position();
01618 if (center.dot(p) < 0.) continue;
01619
01620 float dPhi = phi[sl] - _unused[2][i]->position().phi();
01621 dPhi = fabs(fmod(dPhi, PI2));
01622 if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
01623 if (dPhi < WIDTH[sl] * loadWidth) {
01624 outList.append(_unused[2][i]);
01625 continue;
01626 }
01627 dPhi = phi[sl + 1] - _unused[2][i]->position().phi();
01628 dPhi = fabs(fmod(dPhi, PI2));
01629 if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
01630 if (dPhi < WIDTH[sl] * loadWidth)
01631 outList.append(_unused[2][i]);
01632 }
01633 return outList;
01634 }
01635
01636 int
01637 TConformalFinder::crossPointsInConformal(const AList<TSegment> & inList,
01638 HepPoint3D points[12]) const {
01639
01640 static const float confRadius2[] = {4. / ( 8.3 * 8.3),
01641 4. / (16.9 * 16.9),
01642 4. / (21.7 * 21.7),
01643 4. / (31.3 * 31.3),
01644 4. / (36.1 * 36.1),
01645 4. / (44.1 * 44.1),
01646 4. / (50.5 * 50.5),
01647 4. / (58.5 * 58.5),
01648 4. / (64.9 * 64.9),
01649 4. / (72.9 * 72.9),
01650 4. / (79.3 * 79.3),
01651 4. / (87.4 * 87.4)};
01652
01653
01654 AList<TMLink> forLine;
01655 HepPoint3D center;
01656 unsigned n = inList.length();
01657 for (unsigned i = 0; i < n; i++) {
01658 const HepPoint3D & p = inList[i]->position();
01659 forLine.append(new TMLink(NULL, NULL, p));
01660 center += p;
01661 }
01662 center *= 1. / float(n);
01663
01664
01665 TMLine line(forLine);
01666 int err = line.fit();
01667 if (err) {
01668 #ifdef TRKRECO_DEBUG
01669 std::cout << "crossPoints ... failed due to line fit failure" << std::endl;
01670 #endif
01671 HepAListDeleteAll(forLine);
01672 return -1;
01673 }
01674
01675
01676 for (unsigned i = 0; i < 12; i++) {
01677 HepPoint3D p[2];
01678 float c = - line.a() * line.b();
01679 float d = 1. + line.a() * line.a();
01680 float e = sqrt(confRadius2[i] * d - line.b() * line.b());
01681 p[0].setX((c + e) / d);
01682 p[0].setY(line.a() * p[0].x() + line.b());
01683 p[1].setX((c - e) / d);
01684 p[1].setY(line.a() * p[1].x() + line.b());
01685
01686 float mag0 = (center - p[0]).mag();
01687 float mag1 = (center - p[1]).mag();
01688
01689 if (mag0 < mag1) points[i] = p[0];
01690 else points[i] = p[1];
01691 #ifdef TRKRECO_DEBUG_DETAIL
01692 std::cout << "0,1 = " << p[0] << ", " << p[1] << std::endl;
01693 #endif
01694 }
01695
01696 #ifdef TRKRECO_DEBUG_DETAIL
01697 std::cout << "... center is : " << center << std::endl;
01698 std::cout << "... cross points are : " << std::endl;
01699 for (unsigned i = 0; i < 12; i++) {
01700 std::cout << " " << i << " : " << points[i] << std::endl;
01701 }
01702 #endif
01703 #ifdef TRKRECO_WINDOW
01704 displayStatus("Conf::crossPoints ... line to salvage for conf plane");
01705 _rphiWindow.append(inList, leda_green);
01706 _rphiWindow.oneShot(line, leda_blue);
01707 #endif
01708
01709 HepAListDeleteAll(forLine);
01710 return 0;
01711 }
01712
01713 AList<TSegment>
01714 TConformalFinder::stereoSegments(const TTrack & t) const {
01715 #ifdef TRKRECO_DEBUG
01716 std::cout << "... finding stereo segments" << std::endl;
01717 #endif
01718
01719 const AList<TSegment> & bases = (const AList<TSegment> &) t.segments();
01720 AList<TSegment> seeds;
01721 unsigned n = bases.length();
01722 if (n == 0) return seeds;
01723
01724
01725 TPoint2D points[12];
01726 int err = crossPoints(t, points);
01727
01728
01729
01730 if (err) {
01731 #ifdef TRKRECO_DEBUG
01732 std::cout << "... stereo segment finding failed" << std::endl;
01733 #endif
01734 return seeds;
01735 }
01736
01737
01738
01739 AList<TSegment> list = pickUpSegments(points, float(_stereoLoadWidth), 2);
01740
01741
01742
01743
01744 TPoint2D dir[5];
01745 for (unsigned i = 0; i < n; i++) {
01746 const TSegment & s = * bases[i];
01747 AList<TMLink> ll = s.links();
01748 unsigned sl = ll[0]->wire()->axialStereoLayerId()/4;
01749 dir[sl] = TPoint2D(s.direction());
01750
01751 }
01752
01753
01754 for (unsigned i = 0; i < 5; i++) {
01755 if (dir[i].mag() < .5) {
01756 unsigned j = i;
01757 while ((j < 5) && (dir[j].mag() < .5))
01758 ++j;
01759 if (j > 4) j = 4;
01760 if (dir[j].mag() < .5) {
01761 j = i;
01762 while ((j > 0) && (dir[j].mag() < .5))
01763 --j;
01764 }
01765 dir[i] = dir[j];
01766 }
01767 }
01768
01769 #ifdef TRKRECO_DEBUG_DETAIL
01770 std::cout << " ... direction :" << std::endl;
01771 for (unsigned i = 0; i < 6; i++)
01772 std::cout << " " << i << " : " << dir[i] << std::endl;
01773 std::cout << " ... direction cos :" << std::endl;
01774 #endif
01775 AList<TSegment> badList;
01776 AList<TSegment> toBeChecked[6];
01777 unsigned nL = list.length();
01778 for (unsigned i = 0; i < 6; ++i)
01779 toBeChecked[i].removeAll();
01780
01781 for (unsigned j = 0; j < nL; ++j) {
01782 TSegment & s = * list[j];
01783 AList<TMLink> ll = s.links();
01784 unsigned sl = ll[0]->wire()->axialStereoLayerId()/4;
01785 toBeChecked[sl].append(s);
01786 }
01787
01788 for (unsigned i = 0; i < 6; ++i) {
01789 unsigned nToBeChecked = toBeChecked[i].length();
01790 if (nToBeChecked < 4) continue;
01791 unsigned nBadList = badList.length();
01792 for (unsigned j = 0; j < nToBeChecked; ++j) {
01793 if((badList.length() - nBadList) >= nToBeChecked - 4) break;
01794 TSegment & s = * toBeChecked[i][j];
01795 TPoint2D sDir = s.direction();
01796 unsigned sl = i;
01797 unsigned aSl;
01798
01799 if (sl < 2) aSl = 0;
01800 else if (sl < 4) aSl = 2;
01801 else aSl = 3;
01802 if (dir[aSl].dot(sDir) < 0.7) badList.append(s);
01803 if (sl == 2 || sl == 3 || sl == 4) {
01804 if (dir[aSl - 1].dot(sDir) < 0.6) badList.append(s);
01805 }
01806 else {
01807 if (dir[aSl + 1].dot(sDir) < 0.6) badList.append(s);
01808 }
01809 }
01810 toBeChecked[i].remove(badList);
01811 if (toBeChecked[i].length() > 4) {
01812 AList<TSegment> tmp;
01813 unsigned nSegs = toBeChecked[i].length();
01814 for (unsigned j = 0; j < nSegs; ++j) {
01815 TSegment & s = * toBeChecked[i][j];
01816 if (s.links().length() > 3) tmp.append(s);
01817 }
01818 toBeChecked[i].remove(tmp);
01819 for (unsigned j = 0; j < tmp.length(); ++j) {
01820 AList<TMLink> tmpLinks = tmp[j]->links();
01821 for (unsigned k = 0; k < toBeChecked[i].length(); ++k) {
01822 TSegment & ss = * toBeChecked[i][k];
01823 AList<TMLink> threeLinks = ss.links();
01824 for (unsigned kk = 0; kk < threeLinks.length(); ++kk) {
01825 TMLink & ll = * threeLinks[kk];
01826 if (tmpLinks.hasMember(ll)) {
01827 badList.append(ss);
01828 break;
01829 }
01830 }
01831 }
01832 }
01833 }
01834
01835 toBeChecked[i].remove(badList);
01836 if (toBeChecked[i].length() > 4) {
01837 AList<TSegment> temp;
01838 unsigned nSegs = toBeChecked[i].length();
01839 for(unsigned k=0; k < nSegs; ++k){
01840 TSegment &s = *toBeChecked[i][k];
01841 temp.append(s);
01842 }
01843
01844 for(unsigned y=0; y<nSegs; y++){
01845 for(unsigned x=1; x<nSegs-y; x++){
01846 TSegment &s1 = *temp[x];
01847 TSegment &s2 = *temp[x-1];
01848 if(s1.links().length() < s2.links().length()){
01849 TSegment & s0 = *temp[x-1];
01850 *temp[x-1] = *temp[x];
01851 *temp[x] = s0;
01852 }
01853 }
01854 }
01855 for(unsigned j=4; j<nSegs; ++j){
01856 TSegment &s4 = *temp[j];
01857 badList.append(s4);
01858 }
01859 }
01860
01861
01862
01863 toBeChecked[i].remove(badList);
01864
01865
01866 for(unsigned j = 0; j < toBeChecked[i].length(); ++j){
01867 TSegment &s5 = *toBeChecked[i][j];
01868
01869 }
01870
01871 }
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905 list.remove(badList);
01906
01907
01908 #ifdef TRKRECO_DEBUG
01909 std::cout << " ... bad segments :" << std::endl;
01910 for (unsigned i = 0; i < badList.length(); i++)
01911 badList[i]->dump("hits sort", " ");
01912 #endif
01913
01914 return list;
01915 }
01916
01917 HepVector3D
01918 TConformalFinder::direction(const TSegment & segment) const {
01919 AList<TMLink> forLine;
01920 const TSegment * s = & segment;
01921 while (1) {
01922 if (s->outerLinks().length() != 1) break;
01923 const TSegment * o = s->outerLinks()[0];
01924 const HepPoint3D & p = o->position();
01925 forLine.append(new TMLink(NULL, NULL, p));
01926 s = o;
01927 }
01928
01929 if (forLine.length() == 0)
01930 return segment.direction();
01931 else if (forLine.length() < 2) {
01932 HepAListDeleteAll(forLine);
01933 return segment.direction();
01934 }
01935
01936
01937 TMLine line(forLine);
01938 int err = line.fit();
01939 if (err) {
01940 #ifdef TRKRECO_DEBUG
01941 std::cout << "direction ... failed due to line fit failure" << std::endl;
01942 #endif
01943 HepAListDeleteAll(forLine);
01944 return segment.direction();
01945 }
01946
01947 HepAListDeleteAll(forLine);
01948 Vector3 tmp(-1., -(line.a() + line.b()));
01949 tmp.unit();
01950 return HepVector3D(tmp);
01951 }
01952
01953
01954 void
01955 TConformalFinder::salvage(TTrack & track,
01956 unsigned asSwitch,
01957 const AList<TSegment> & bads) const {
01958 #ifdef TRKRECO_DEBUG_DETAIL
01959 std::cout << "... salvaging in real plane" << std::endl;
01960 #endif
01961
01962
01963 TPoint2D points[12];
01964 int err = crossPoints(track, points);
01965 if (err) {
01966 #ifdef TRKRECO_DEBUG_DETAIL
01967 std::cout << "... salvaging failed in calculation of intersection"
01968 << std::endl;
01969 #endif
01970 return;
01971 }
01972
01973
01974 AList<TSegment> toBeChecked = pickUpSegments(points,
01975 float(_salvageLoadWidth),
01976 asSwitch);
01977 toBeChecked.remove(bads);
01978 toBeChecked.remove(track.segments());
01979 toBeChecked = trackSide(track, toBeChecked);
01980
01981
01982 AList<TMLink> links;
01983 unsigned nB = toBeChecked.length();
01984 for (unsigned i = 0; i < nB; i++) {
01985 const AList<TMLink> & tmp = toBeChecked[i]->links();
01986 unsigned nL = tmp.length();
01987 for (unsigned j = 0; j < nL; j++) {
01988 if (tmp[j]->hit()->track()) continue;
01989 links.append(tmp[j]);
01990 }
01991 }
01992
01993
01994 AList<TMLink> all = pickUpLinks(points, float(_salvageLoadWidth), asSwitch);
01995 #ifdef TRKRECO_DEBUG
01996 std::cout<<" pickUpLinks length = "<<all.length()<<std::endl;
01997 #endif
01998
01999 all = trackSide(track, all);
02000 all.remove(track.links());
02001 #ifdef TRKRECO_DEBUG
02002 std::cout<<" after reomove pickUpLinks length = "<<all.length()<<std::endl;
02003 #endif
02004 links.append(all);
02005
02006 #ifdef TRKRECO_WINDOW
02007 AList<TMLink> tmp = links;
02008 #endif
02009
02010
02011 _builder.salvage(track, links);
02012
02013
02014 const AList<TMLink> & usedLinks = track.links();
02015 for (unsigned i = 0; i < nB; i++) {
02016 TSegment & segment = * toBeChecked[i];
02017 AList<TMLink> used = Links(segment, track);
02018 if (used.length() == 0) continue;
02019
02020 track.segments().append(segment);
02021 segment.tracks().append(track);
02022 }
02023
02024 #ifdef TRKRECO_WINDOW
02025 displayStatus("Conf::salvage ... results");
02026 TTrackBase y(tmp);
02027 _rphiWindow.append(y, leda_red);
02028 _rphiWindow.append(toBeChecked, leda_red);
02029 _rphiWindow.oneShot(track, leda_green);
02030 #endif
02031
02032
02033 return;
02034 }
02035
02036 int
02037 TConformalFinder::crossPoints(const TTrack & track,
02038 TPoint2D points[12]) const {
02039
02040 #ifdef TRKRECO_DEBUG_DETAIL
02041 std::cout << " cross points are :" << std::endl;
02042 #endif
02043
02044
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070
02071 static const float R[] = { 7.3,
02072 12.21,
02073 18.9,
02074 25.375,
02075 31.86,
02076 39.16,
02077 45.685,
02078 52.215,
02079 58.665,
02080 65.91,
02081 72.37,
02082 80.};
02083 static const float R2[] = { R[0]*R[0],
02084 R[1]*R[1],
02085 R[2]*R[2],
02086 R[3]*R[3],
02087 R[4]*R[4],
02088 R[5]*R[5],
02089 R[6]*R[6],
02090 R[7]*R[7],
02091 R[8]*R[8],
02092 R[9]*R[9],
02093 R[10]*R[10],
02094 R[11]*R[11]};
02095
02096
02097 static const TPoint2D o(0., 0.);
02098 TPoint2D c = track.center();
02099 TPoint2D co = - c;
02100 double r = fabs(track.radius());
02101 double r2 = r * r;
02102 double d2 = c.mag2();
02103 double d = sqrt(d2);
02104 double sl = - c.x() / c.y();
02105
02106 unsigned nOk = 0;
02107 for (unsigned i = 0; i < 12; i++) {
02108 double minR = r < R[i] ? r : R[i];
02109 double maxR = r < R[i] ? R[i] : r;
02110
02111 #ifdef TRKRECO_DEBUG_DETAIL
02112 std::cout << " r,R ... " << minR << ", " << maxR << ", " << d << std::endl;
02113 #endif
02114
02115 if ((r + R[i] < d) || (minR + d < maxR)) {
02116 points[i] = o;
02117 continue;
02118 }
02119 ++nOk;
02120 double a = R2[i] + d2 - r2;
02121 double s = sqrt(4. * R2[i] * d2 - a * a);
02122 double q = 0.5 * a / c.y();
02123 points[i].x(0.5 * (c.x() * a + c.y() * s) / d2);
02124 points[i].y(q + sl * points[i].x());
02125 const double Bz = -1000*(m_pmgnIMF->getReferField());
02126 if (co.cross(points[i] - c) * track.charge() * Bz > 0.) {
02127 points[i].x(0.5 * (c.x() * a - c.y() * s) / d2);
02128 points[i].y(q + sl * points[i].x());
02129 }
02130 #ifdef TRKRECO_DEBUG_DETAIL
02131 std::cout << " " << i << " : " << points[i] << std::endl;
02132 std::cout << " chg=" << track.charge()<<std::endl;
02133 std::cout << ", c=" << c << ", co=" <<std::endl;
02134 std::cout << ", " << co.cross(points[i] - c) * track.charge() << std::endl;
02135 std::cout << " " << 0.5 * (c.x() * a + c.y() * s) / d2;
02136 std::cout << ", " << q + sl * (0.5 * (c.x() * a + c.y() * s) / d2);
02137 std::cout << " " << 0.5 * (c.x() * a - c.y() * s) / d2;
02138 std::cout << ", " << q + sl * (0.5 * (c.x() * a - c.y() * s) / d2);
02139 std::cout << std::endl;
02140 #endif
02141
02142 }
02143 if (nOk) return 0;
02144 else return -1;
02145 }
02146
02147 AList<TSegment>
02148 TConformalFinder::pickUpSegments(const TPoint2D x[12],
02149 float loadWidth,
02150 unsigned axialStereoSwitch) const {
02151 static const TPoint2D O(0., 0.);
02152 AList<TSegment> outList;
02153
02154
02155
02156 for (unsigned sl = 0; sl < 2; sl++) {
02157 unsigned nMax = 5;
02158 if (sl == 0) {
02159 if (! (axialStereoSwitch & 1)) continue;
02160 }
02161 else {
02162 nMax = 6;
02163 if (! (axialStereoSwitch & 2)) continue;
02164 }
02165 for (unsigned i = 0; i < nMax; i++) {
02166
02167 unsigned layerId;
02168 switch(sl){
02169 case 0:
02170 switch(i){
02171 case 0:layerId = 2;break;
02172 case 1:layerId = 3;break;
02173 case 2:layerId = 4;break;
02174 case 3:layerId = 9;break;
02175 case 4:layerId = 10;break;
02176 }
02177 case 1:
02178 switch(i){
02179 case 0:layerId = 0;break;
02180 case 1:layerId = 1;break;
02181 case 2:layerId = 5;break;
02182 case 3:layerId = 6;break;
02183 case 4:layerId = 7;break;
02184 case 5:layerId = 8;break;
02185 }
02186 }
02187
02188
02189 if (x[layerId] == O) continue;
02190 float a = WIDTH[layerId] * loadWidth;
02191 float phi0 = x[layerId].phi();
02192 float phi1 = x[layerId+1].phi();
02193 float phiC = (phi0 + phi1) / 2.;
02194 float phiD = fabs(phi0 - phiC);
02195 if (phiD > M_PI_2) {
02196 phiC += M_PI;
02197 phiD = M_PI - phiD;
02198 }
02199 phiD += a;
02200
02201 const AList<TSegment> & list = _allUnused[sl][i];
02202 unsigned n = list.length();
02203 #ifdef TRKRECO_DEBUG_DETAIL
02204 std::cout << " pickUpSegments ... phi0,1,C,D=" << phi0 << ",";
02205 std::cout << phi1 << "," << phiC << "," << phiD <<" nhits "<<n<< std::endl;
02206 #endif
02207 for (unsigned j = 0; j < n; j++) {
02208 TSegment & s = * list[j];
02209
02210 #ifdef TRKRECO_DEBUG_DETAIL
02211 s.dump("hits sort", " ");
02212 #endif
02213
02214 bool found = false;
02215 const AList<TMLink> & inners = s.inners();
02216 unsigned m = inners.length();
02217 for (unsigned k = 0; k < m; k++) {
02218 float phi = inners[k]->position().phi();
02219 if (phi < 0.) phi += PI2;
02220 float dPhi = fabs(phi - phiC);
02221 if (dPhi > M_PI) dPhi = PI2 - dPhi;
02222 #ifdef TRKRECO_DEBUG_DETAIL
02223 std::cout << " " << inners[k]->wire()->name();
02224 std::cout << ", phi,dPhi=" << phi << "," << dPhi << std::endl;
02225 #endif
02226 if (dPhi < phiD) {
02227 outList.append(s);
02228 found = true;
02229 break;
02230 }
02231 }
02232
02233 if (found) continue;
02234 const AList<TMLink> & outers = s.outers();
02235 m = outers.length();
02236 for (unsigned k = 0; k < m; k++) {
02237 float phi = outers[k]->position().phi();
02238 if (phi < 0.) phi += PI2;
02239 float dPhi = fabs(phi - phiC);
02240 if (dPhi > M_PI) dPhi = PI2 - dPhi;
02241 #ifdef TRKRECO_DEBUG_DETAIL
02242 std::cout << " " << outers[k]->wire()->name();
02243 std::cout << ", phi,dPhi=" << phi << "," << dPhi << std::endl;
02244 #endif
02245 if (dPhi < phiD) {
02246 outList.append(s);
02247 found = true;
02248 break;
02249 }
02250 }
02251 }
02252 }
02253 }
02254 return outList;
02255 }
02256
02257 AList<TMLink>
02258 TConformalFinder::pickUpLinks(const TPoint2D x[12],
02259 float loadWidth,
02260 unsigned axialStereoSwitch) const {
02261
02262 static const TPoint2D O(0., 0.);
02263 AList<TMLink> outList;
02264
02265 unsigned nBad = _unused[2].length();
02266 for (unsigned i = 0; i < nBad; i++) {
02267 unsigned sl = _unused[2][i]->wire()->superLayerId();
02268 unsigned as = 1;
02269 if (_unused[2][i]->wire()->axial()) as = 0;
02270 if (as == 0) {
02271 if (! (axialStereoSwitch & 1)) continue;
02272 }
02273 else {
02274 if (! (axialStereoSwitch & 2)) continue;
02275 }
02276
02277 float a = WIDTH[sl] * loadWidth;
02278 float phi0 = x[sl].phi();
02279 float phi1 = x[sl + 1].phi();
02280 float phi = _unused[2][i]->position().phi();
02281
02282 if (phi < 0.) phi += PI2;
02283 if (phi1 < phi0) {
02284 phi1 = phi0;
02285 phi0 = x[sl + 1].phi();
02286 }
02287 float dPhi = phi1 - phi0;
02288 if (dPhi < M_PI) {
02289 phi0 -= a;
02290 phi1 += a;
02291 if (phi > phi0 && phi < phi1) outList.append(_unused[2][i]);
02292 }
02293 else {
02294 phi0 += a;
02295 phi1 -= a;
02296 if (phi < phi0 || phi > phi1) outList.append(_unused[2][i]);
02297 }
02298 #ifdef TRKRECO_DEBUG_DETAIL
02299 std::cout << "TConformalFinder:: pickUpLinks " << std::endl;
02300 std::cout << "a "<<a <<" phi0 "<<phi0<<" phi1 "<<phi1<<" phi "<< phi<<" dPhi "<<dPhi << std::endl;
02301 #endif
02302 }
02303 return outList;
02304 }
02305
02306 void
02307 TConformalFinder::slowFinding2D(unsigned level) {
02308
02309 #ifdef TRKRECO_DEBUG_DETAIL
02310 std::cout << "... slow finding (2D) : level = " << level << std::endl;
02311 #endif
02312 #ifdef TRKRECO_DEBUG
02313 _stage = ConformalSlow2D;
02314 #endif
02315
02316 AList<TSegment> * segments = & _allUnused[0][0];
02317
02318 unsigned id = _2DTracks.length() + _3DTracks.length();
02319 unsigned seedLayer = 5;
02320 while (seedLayer--) {
02321
02322
02323 AList<TSegment> tmp = segments[seedLayer];
02324 unsigned n = tmp.length();
02325 for (unsigned i = 0; i < n; i++) {
02326 AList<TTrack> trackCandidates;
02327 TSegment & current = * tmp[i];
02328 if (current.links().length() < 3) continue;
02329 if (current.innerLinks().length() != 1) continue;
02330 if (current.innerLinks()[0]->links().length() < 3) continue;
02331 AList<TSegment> seeds;
02332 seeds.append(current);
02333 seeds.append(current.innerLinks()[0]);
02334
02335 std::string name = "s" + itostring(id);
02336
02337 #ifdef TRKRECO_DEBUG
02338 std::cout << "==> slow building : " << name << std::endl;
02339 #endif
02340
02341
02342 for (int k = 0; k < seeds.length(); ++k)
02343 seeds[k]->solveLR();
02344
02345 TTrack * t = expand(seeds);
02346 if (t) {
02347 AList<TSegment> bads;
02348 t->name(name);
02349 salvage(* t, 1, bads);
02350 if (! trackQuality(* t)) {
02351 deleteTrack(* t);
02352 continue;
02353 }
02354 t->finder(TrackSlowFinder);
02355 t->quality(TrackQuality2D);
02356 t->fitting(TrackFitGlobal);
02357 trackCandidates.append(t);
02358
02359 #ifdef TRKRECO_DEBUG
02360 std::cout << "... 2D finished : " << t->name() << std::endl;
02361 t->dump("hits sort flag pull", " ");
02362 #endif
02363 #ifdef TRKRECO_WINDOW
02364 displayStatus("Conf::expand ... finished");
02365 _rphiWindow.oneShot(* t, leda_green);
02366 #endif
02367 removeUsedSegments(trackCandidates);
02368 _2DTracks.append(trackCandidates);
02369 id = _2DTracks.length() + _3DTracks.length();
02370 }
02371 }
02372
02373
02374
02375
02376
02377 }
02378 }
02379
02380 TTrack *
02381 TConformalFinder::expand(AList<TSegment> & seeds) const {
02382 #ifdef TRKRECO_WINDOW
02383 displayStatus("Conf::expand ... seeds");
02384 _rphiWindow.oneShot(seeds, leda_green);
02385 #endif
02386 #ifdef TRKRECO_DEBUG
02387 std::cout << "... expand : # of seeds = " << seeds.length() << std::endl;
02388 #endif
02389 TTrack * t = NULL;
02390 if (seeds.length() > 2) {
02391 TTrack * t = _builder.buildRphi(seeds);
02392 if (t) {
02393 if (! trackQuality(* t)) {
02394 deleteTrack(* t);
02395 t = NULL;
02396 }
02397 }
02398 }
02399 if (t == NULL) {
02400 AList<TMLink> links = Links(seeds);
02401 TCircle c(links);
02402 c.fit();
02403 t = new TTrack(c);
02404 t->fit();
02405 t->segments().append(seeds);
02406 }
02407 return expand(* t);
02408 }
02409
02410 TTrack *
02411 TConformalFinder::expand(TTrack & ti) const {
02412 #ifdef TRKRECO_WINDOW
02413 displayStatus("Conf::expand ... seed track");
02414 _rphiWindow.oneShot(ti, leda_green);
02415 #endif
02416 #ifdef TRKRECO_DEBUG
02417 std::cout << "... expand : by a track" << std::endl;
02418 #endif
02419
02420 TTrack * t = & ti;
02421 AList<TSegment> seeds = t->segments();
02422 unsigned nSegments = seeds.length();
02423 unsigned nCores = t->nCores();
02424 unsigned nTrial = 0;
02425 while (1) {
02426
02427
02428 unsigned sl = SuperLayer(t->cores());
02429 unsigned inner;
02430 unsigned outer;
02431 targetSuperLayer(sl, inner, outer);
02432
02433 unsigned target = inner;
02434 if (target == 11)
02435 target = outer;
02436 if (target == 11)
02437 break;
02438
02439
02440 AList<TSegment> segments = targetSegments(* t, target);
02441 if (segments.length() == 0) {
02442 unsigned fl[5] = {2, 3, 4, 9, 10};
02443 unsigned inner2 = 0;
02444 unsigned outer2 = 5;
02445 for (int i = 0; i < 5; i++) {
02446 if (fl[i] == inner) inner2 = i;
02447 if (fl[i] == outer) outer2 = i;
02448 }
02449 if (inner > 2&& inner2>0) target = fl[inner2 - 1];
02450 if (target == 9) {
02451 target = fl[outer2 + 1];
02452 }
02453 if (target == 13) break;
02454 segments = targetSegments(* t, target);
02455 }
02456 #ifdef TRKRECO_DEBUG
02457 std::cout << "... expand : segments to be checked = ";
02458 std::cout << segments.length() << std::endl;
02459 #endif
02460 #ifdef TRKRECO_WINDOW
02461 displayStatus("Conf::expand ... candidate segments");
02462 _rphiWindow.append(segments, leda_blue);
02463 _rphiWindow.oneShot(* t, leda_green);
02464 #endif
02465
02466
02467 unsigned n = segments.length();
02468 AList<TTrack> tracks;
02469 for (unsigned i = 0; i < n; i++) {
02470 segments[i]->solveLR();
02471 AList<TSegment> seed0 = seeds;
02472 seed0.append(segments[i]);
02473 TTrack * t0 = _builder.buildRphi(seed0);
02474 if (! t0) continue;
02475 if ((t0->segments().length() <= nSegments) ||
02476 (t0->nCores() <= nCores) ||
02477 (SuperLayer(t0->cores()) <= sl)) {
02478 deleteTrack(* t0);
02479 continue;
02480 #ifdef TRKRECO_DEBUG
02481 std::cout << "... expand : candidate deleted" << std::endl;
02482 #endif
02483 }
02484 tracks.append(t0);
02485 }
02486
02487 #ifdef TRKRECO_DEBUG
02488 std::cout << "... expand : target superlayer = " << target;
02489 std::cout << ", # of track candidates = " << tracks.length()
02490 << std::endl;
02491 #endif
02492
02493
02494 if (tracks.length() == 0) {
02495 AList<TMLink> links = targetLinks(* t, target);
02496 _builder.salvage(* t, links);
02497 unsigned newSl = SuperLayer(t->cores());
02498 if (newSl & (1 << (target))) {
02499 #ifdef TRKRECO_DEBUG
02500 std::cout << "... expand : links to be appended = ";
02501 std::cout << links.length() << std::endl;
02502 #endif
02503 #ifdef TRKRECO_WINDOW
02504 displayStatus("Conf::expand ... candidate links");
02505 _rphiWindow.append(links, leda_blue);
02506 _rphiWindow.oneShot(* t, leda_green);
02507 #endif
02508 TTrack * t0 = new TTrack(* t);
02509 if (! t0) break;
02510 if ((t0->nCores() <= nCores) ||
02511 (SuperLayer(t0->cores()) <= sl)) {
02512 deleteTrack(* t0);
02513 break;
02514 #ifdef TRKRECO_DEBUG
02515 std::cout << "... expand : candidate deleted" << std::endl;
02516 #endif
02517 }
02518 tracks.append(t0);
02519 }
02520 else
02521 break;
02522 }
02523
02524 #ifdef TRKRECO_DEBUG
02525 std::cout << "... expand : # of track candidates = " << tracks.length()
02526 << std::endl;
02527 #endif
02528
02529
02530 unsigned nTracks = tracks.length();
02531 if (nTracks == 0) break;
02532 TTrack * tNew = tracks[0];
02533 unsigned nBestCores = tNew->nCores();
02534 for (unsigned i = 1; i < nTracks; i++) {
02535 if (tracks[i]->nCores() > nBestCores) {
02536 deleteTrack(* tNew);
02537 tNew = tracks[i];
02538 nBestCores = tNew->nCores();
02539 }
02540 else {
02541 deleteTrack(* tracks[i]);
02542 }
02543 }
02544 nSegments = tNew->segments().length();
02545 nCores = nBestCores;
02546 seeds = tNew->segments();
02547
02548
02549 if (! trackQuality(* tNew)) {
02550 deleteTrack(* tNew);
02551 break;
02552 }
02553 deleteTrack(* t);
02554 t = tNew;
02555 }
02556
02557 #ifdef TRKRECO_DEBUG
02558 if (t == & ti) std::cout << "... expand : failed" << std::endl;
02559 else std::cout << "... expand : track expanded" << std::endl;
02560 #endif
02561 #ifdef TRKRECO_WINDOW
02562 displayStatus("Conf::expand ... finished");
02563 _rphiWindow.oneShot(* t, leda_green);
02564 #endif
02565
02566 return t;
02567 }
02568
02569 void
02570 TConformalFinder::targetSuperLayer(unsigned sl,
02571 unsigned & inner,
02572 unsigned & outer) const {
02573 inner = 11;
02574 outer = 11;
02575 bool innerFound = false;
02576 bool outerFound = false;
02577 for (unsigned i = 0; i < 5; i++) {
02578 unsigned fl[5] ={2, 3, 4, 9, 10};
02579 if (! innerFound) {
02580 if (sl & (1 << (fl[i]))) innerFound = true;
02581 else inner = fl[i];
02582 }
02583 if (! outerFound) {
02584 if (sl & (1 << (fl[4 - i]))) outerFound = true;
02585 else outer = fl[4 - i];
02586 }
02587 }
02588 }
02589
02590 AList<TSegment>
02591 TConformalFinder::targetSegments(const TTrack & t, unsigned sl) const {
02592 AList<TSegment> targets;
02593
02594
02595 TPoint2D points[12];
02596 int err = crossPoints(t, points);
02597 if (err) {
02598 #ifdef TRKRECO_DEBUG
02599 std::cout << " ... no target found" << std::endl;
02600 #endif
02601 return targets;
02602 }
02603
02604
02605 AList<TSegment> toBeChecked = pickUpSegments(points, 3, 1);
02606 unsigned n = toBeChecked.length();
02607 for (unsigned i = 0; i < n; i++) {
02608 if (toBeChecked[i]->superLayerId() == sl)
02609 targets.append(toBeChecked[i]);
02610 }
02611
02612 return targets;
02613 }
02614
02615 AList<TMLink>
02616 TConformalFinder::targetLinks(const TTrack & t, unsigned sl) const {
02617 AList<TMLink> targets;
02618
02619
02620 TPoint2D points[12];
02621 int err = crossPoints(t, points);
02622 if (err) {
02623 #ifdef TRKRECO_DEBUG
02624 std::cout << " ... no target found" << std::endl;
02625 #endif
02626 return targets;
02627 }
02628
02629
02630 AList<TMLink> toBeChecked = pickUpLinks(points, 3, 1);
02631 unsigned n = toBeChecked.length();
02632 for (unsigned i = 0; i < n; i++) {
02633 if (toBeChecked[i]->wire()->superLayerId() == sl)
02634 targets.append(toBeChecked[i]);
02635 }
02636
02637 return targets;
02638 }
02639
02640 AList<TSegment>
02641 TConformalFinder::refineSegments(const TTrack & t) const {
02642 const AList<TSegment> & original = t.segments();
02643 AList<TSegment> outList;
02644 unsigned n = original.length();
02645 for (unsigned i = 0; i < n; i++) {
02646 AList<TMLink> tmp = Links(* original[i], t);
02647 if (tmp.length() >= _minNLinksForSegmentInRefine)
02648 outList.append(original[i]);
02649 }
02650 unsigned sl = SuperLayer(outList);
02651 unsigned nCMax = 0;
02652 unsigned nStart = 0;
02653 unsigned nC = 0;
02654 unsigned nS = 0;
02655 unsigned fl[5] = {2, 3, 4, 9, 10};
02656 for (unsigned i = 0; i < 5; i++) {
02657 if (sl & (1 << fl[i])) {
02658 ++nC;
02659 if (nC == 1) nS = i;
02660 if (nC > nCMax) {
02661 nCMax = nC;
02662 nStart = nS;
02663 }
02664 }
02665 else {
02666 nC = 0;
02667 }
02668 }
02669
02670
02671
02672 outList.removeAll();
02673 if (nCMax >= _minNSegments) {
02674 for (unsigned i = 0; i < n; i++) {
02675 unsigned id = original[i]->superLayerId();
02676 if ((id >= fl[nStart]) && (id < fl[nStart + nCMax]))
02677 outList.append(original[i]);
02678 }
02679 }
02680
02681 #ifdef TRKRECO_DEBUG
02682 std::cout << "... refine segments : orignal sl = ";
02683 bitDisplay(sl, 11, 0);
02684 std::cout << ", output sl = ";
02685 bitDisplay(SuperLayer(outList), 11, 0);
02686 std::cout << std::endl;
02687 #endif
02688
02689 return outList;
02690 }
02691
02692 bool
02693 TConformalFinder::trackQuality(const TTrack & t) const {
02694 unsigned n1 = NSuperLayers(t.links());
02695 unsigned n3 = NSuperLayers(t.links(), 2);
02696
02697
02698 if(fabs(t.radius())<15.) return false;
02699
02700 #ifdef TRKRECO_WINDOW
02701 std::cout << "... trackQuality : n1,n3,nMissing=" << n1 << "," << n3;
02702 std::cout << "," << NMissingAxialSuperLayers(t.links()) << std::endl;
02703 #endif
02704 #ifdef TRKRECO_WINDOW
02705 displayStatus("trackQuality");
02706 if ((n3 < 2) || (NMissingAxialSuperLayers(t.links()) > 1))
02707 _rphiWindow.oneShot(t, leda_red);
02708 else
02709 _rphiWindow.oneShot(t, leda_green);
02710 #endif
02711 if (n3 < 2) return false;
02712
02713
02714 return true;
02715 }
02716
02717 void
02718 TConformalFinder::refineLinks(TTrack & t, unsigned minN) const {
02719 const AList<TMLink> & links = t.links();
02720 AList<TMLink> sl[11];
02721 unsigned n = links.length();
02722 for (unsigned i = 0; i < n; i++)
02723 sl[links[i]->wire()->superLayerId()].append(links[i]);
02724 AList<TMLink> toBeRemoved;
02725 for (unsigned i = 0; i < 11; i++)
02726 if (sl[i].length() < minN)
02727 toBeRemoved.append(sl[i]);
02728 t.remove(toBeRemoved);
02729 #ifdef TRKRECO_DEBUG
02730 std::cout << "... refining links : removed links = " << toBeRemoved.length();
02731 std::cout << std::endl;
02732 #endif
02733 }
02734
02735 AList<TMLink>
02736 TConformalFinder::trackSide(const TTrack & t, const AList<TMLink> & a) const {
02737 static const TPoint2D o(0., 0.);
02738 TPoint2D c = t.center();
02739 TPoint2D co = - c;
02740
02741 AList<TMLink> tmp;
02742 unsigned n = a.length();
02743 for (unsigned i = 0; i < n; i++) {
02744 TPoint2D x = a[i]->wire()->xyPosition();
02745 const double Bz = -1000*(m_pmgnIMF->getReferField());
02746 if (co.cross(x - c) * t.charge() * Bz < 0.)
02747 tmp.append(a[i]);
02748 }
02749 return tmp;
02750 }
02751
02752 AList<TSegment>
02753 TConformalFinder::trackSide(const TTrack & t,
02754 const AList<TSegment> & a) const {
02755 static const TPoint2D o(0., 0.);
02756 TPoint2D c = t.center();
02757 TPoint2D co = - c;
02758
02759 AList<TSegment> tmp;
02760 unsigned n = a.length();
02761 for (unsigned i = 0; i < n; i++) {
02762 const AList<TMLink> & b = a[i]->links();
02763 bool bad = false;
02764 unsigned nB = b.length();
02765 for (unsigned j = 0; j < nB; j++) {
02766 TPoint2D x = b[j]->wire()->xyPosition();
02767
02768 const double Bz = -1000*(m_pmgnIMF->getReferField());
02769 if (co.cross(x - c) * t.charge()* Bz > 0.) {
02770 bad = true;
02771 break;
02772 }
02773 }
02774 if (bad) continue;
02775 tmp.append(a[i]);
02776 }
02777 return tmp;
02778 }
02779
02780 void
02781 TConformalFinder::resolveHits(AList<TTrack> & tracks) const {
02782 unsigned nTracks = tracks.length();
02783 if (nTracks < 2) return;
02784
02785 for (unsigned i = 0; i < nTracks - 1; i++) {
02786 TTrack & t0 = * tracks[i];
02787 const AList<TMLink> & links0 = t0.links();
02788 unsigned n0 = links0.length();
02789 AList<TMLink> remove0;
02790 for (unsigned j = i + 1; j < nTracks; j++) {
02791 TTrack & t1 = * tracks[j];
02792 const AList<TMLink> & links1 = t1.links();
02793 unsigned n1 = links1.length();
02794 AList<TMLink> remove1;
02795
02796
02797 for (unsigned k = 0; k < n0; k++) {
02798 TMLink & l = * links0[k];
02799
02800
02801 if (links1.hasMember(l)) {
02802 TMLink l0(NULL, l.hit());
02803 TMLink l1(NULL, l.hit());
02804 t0.approach(l0);
02805 t1.approach(l1);
02806 if (l0.pull() > l1.pull())
02807 remove0.append(l);
02808 else
02809 remove1.append(l);
02810 }
02811 }
02812
02813 if (remove1.length()) {
02814 t1.remove(remove1);
02815 t1.fit();
02816 }
02817 }
02818
02819 if (remove0.length()) {
02820 t0.remove(remove0);
02821 t0.fit();
02822 }
02823 }
02824 }
02825
02826 AList<TSegment>
02827 TConformalFinder::stereoSegmentsFromBadHits(const TTrack & t) const {
02828
02829 AList<TSegment> output;
02830 for (unsigned i = 0; i < 6; i++)
02831 output.append(new TSegment());
02832
02833
02834 TPoint2D points[12];
02835 int err = crossPoints(t, points);
02836 if (err) {
02837 #ifdef TRKRECO_DEBUG_DETAIL
02838 std::cout << "... conf::stereoSegBadHits : "
02839 << "error in calculation of intersection" << std::endl;
02840 #endif
02841 return output;
02842 }
02843
02844 static const TPoint2D O(0., 0.);
02845 unsigned nBad = _unused[2].length();
02846 for (unsigned i = 0; i < nBad; i++) {
02847 unsigned sl = _unused[2][i]->wire()->superLayerId();
02848 unsigned asl = _unused[2][i]->wire()->axialStereoLayerId()/4;
02849 unsigned as;
02850 if (_unused[2][i]->wire()->axial())
02851 as = 0;
02852 else as = 1;
02853
02854 if (as == 0) continue;
02855
02856 float a = WIDTH[sl] * _salvageLoadWidth;
02857 float phi0 = points[sl].phi();
02858 float phi1 = points[sl + 1].phi();
02859 float phi = _unused[2][i]->position().phi();
02860 if (phi < 0.) phi += PI2;
02861 if (phi1 < phi0) {
02862 phi1 = phi0;
02863 phi0 = points[sl + 1].phi();
02864 }
02865 float dPhi = phi1 - phi0;
02866 if (dPhi < M_PI) {
02867 phi0 -= a;
02868 phi1 += a;
02869 if (phi > phi0 && phi < phi1)
02870
02871 output[asl / 4]->append(* _unused[2][i]);
02872 }
02873 else {
02874 phi0 += a;
02875 phi1 -= a;
02876 if (phi < phi0 || phi > phi1)
02877 output[asl / 4]->append(* _unused[2][i]);
02878 }
02879 }
02880
02881 return output;
02882 }
02883
02884 void
02885 TConformalFinder::findSegmentsPerfect(void) {
02886
02887
02888 AList<TMLink> links[11];
02889 unsigned nHits = _hits[2].length();
02890 for (unsigned i = 0; i < nHits; i++) {
02891 TMLink & l = * _hits[2][i];
02892 links[l.wire()->superLayerId()].append(l);
02893 }
02894
02895
02896 const unsigned nHep = TTrackHEP::list().length();
02897
02898
02899 for (unsigned i = 0; i < 11; i++) {
02900 unsigned superLayerId = i / 2;
02901 unsigned id = i % 2;
02902
02903
02904 AList<TMLink> ** seeds =
02905 (AList<TMLink> **) malloc(nHep * sizeof(AList<TMLink> *));
02906 for (unsigned j = 0; j < nHep; j++)
02907 seeds[j] = NULL;
02908
02909
02910 unsigned n = links[i].length();
02911 for (unsigned j = 0; j < n; j++) {
02912 TMLink & l = * links[i][j];
02913 const TMDCWireHitMC * mc = l.hit()->mc();
02914 if (! l.hit()->mc()) {
02915 std::cout << "TConformalFinder::findSegmentsPerfect !!! "
02916 << "no MC info. found ... aborted" << std::endl;
02917 return;
02918 }
02919 const unsigned id = mc->hep()->id();
02920
02921 if (! seeds[id])
02922 seeds[id] = new AList<TMLink>();
02923 seeds[id]->append(l);
02924 }
02925
02926
02927 AList<TSegment> segments;
02928 for (unsigned j = 0; j < nHep; j++) {
02929 if (! seeds[j]) continue;
02930
02931 const unsigned length = seeds[j]->length();
02932 if (length < _minNLinksForSegment) continue;
02933 if (length > _maxNLinksForSegment) continue;
02934 if (Width(* seeds[j]) > _maxWidthForSegment) continue;
02935
02936 TSegment & s = * new TSegment();
02937 s.append(* seeds[j]);
02938 segments.append(s);
02939 }
02940 _allSegments[id][superLayerId].append(segments);
02941 _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
02942
02943
02944 for (unsigned j = 0; j < nHep; j++) {
02945 if (seeds[j])
02946 delete seeds[j];
02947 }
02948 free(seeds);
02949 }
02950
02951 #ifdef TRKRECO_DEBUG_DETAIL
02952 std::cout << "... segment finding(perfect) finished" << std::endl;
02953 #endif
02954 }
02955
02956 void
02957 TConformalFinder::findSegmentsTsf(void) {
02958
02959 AList<TMLink> links[11];
02960 unsigned n = _hits[2].length();
02961 for (unsigned i = 0; i < n; i++) {
02962 TMLink & l = * _hits[2][i];
02963 links[l.wire()->superLayerId()].append(l);
02964 }
02965
02966
02967 TMDCTsf * tsf[11];
02968
02969
02970 for (unsigned i = 0; i < 11; ++i)
02971 tsf[i] = new TMDCTsf(i);
02972
02973 for (unsigned i = 0; i < 11; i++) {
02974 unsigned superLayerId;
02975 unsigned id;
02976 switch(i){
02977 case 0: superLayerId = 0; id = 1; break;
02978 case 1: superLayerId = 1; id = 1; break;
02979 case 2: superLayerId = 0; id = 0; break;
02980 case 3: superLayerId = 1; id = 0; break;
02981 case 4: superLayerId = 2; id = 0; break;
02982 case 5: superLayerId = 2; id = 1; break;
02983 case 6: superLayerId = 3; id = 1; break;
02984 case 7: superLayerId = 4; id = 1; break;
02985 case 8: superLayerId = 5; id = 1; break;
02986 case 9: superLayerId = 3; id = 0; break;
02987 case 10: superLayerId = 4; id = 0; break;
02988 default: break;
02989 }
02990
02991 AList<TSegment> tmp = tsf[i]->segments(links[i]);
02992
02993
02994 _allSegments[id][superLayerId].append(tmp);
02995 _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
02996 delete tsf[i];
02997
02998
02999
03000
03001
03002
03003
03004 }
03005 }