BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
TConformalFinder.cxx
Go to the documentation of this file.
1//-----------------------------------------------------------------------------
2// $Id: TConformalFinder.cxx,v 1.49 2022/01/28 22:04:42 maqm Exp $
3//-----------------------------------------------------------------------------
4// Filename : TConformalFinder.cc
5// Section : Tracking
6// Owner : Yoshi Iwasaki
7// Email : [email protected]
8//-----------------------------------------------------------------------------
9// Description : A class to find tracks with the conformal method.
10// See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
11//-----------------------------------------------------------------------------
12// TConformalLink removed. TSegment added.
13//-----------------------------------------------------------------------------
14
15#include <algorithm>
16#include <string.h>
17#include "CLHEP/String/Strings.h"
18//#include "basf/basfshm.h"
19#include "TrkReco/TMDC.h"
20#include "TrkReco/TMDCUtil.h"
21#include "TrkReco/TTrackManager.h"
22#include "TrkReco/TConformalFinder.h"
23#include "TrkReco/TConformalFinder0.h"
24#include "TrkReco/THistogram.h"
25#include "TrkReco/TMDCTsf.h"
26#include "TrkReco/TCircle.h"
27#include "TrkReco/TTrack.h"
28#include "TrkReco/TMLine.h"
29#include "TrkReco/TTrackHEP.h"
30#include "TrkReco/TMDCWireHitMC.h"
31
32#ifdef TRKRECO_DEBUG
33unsigned TConformalFinder::_stage = ConformalOutside;
34#endif
35#ifdef TRKRECO_WINDOW
36#include "TrkReco/TWindow.h"
37#endif
38
39//extern BasfSharedMem * BASF_Sharedmem;
40
41static const float PI2 = 2. * M_PI;
42static const float WIDTH[] = {PI2 / 40,
43 PI2 / 64, //the first two are stereo, tighter would be better for noise? Liu
44 PI2 / 76,
45 PI2 / 100,
46 PI2 / 128,
47 PI2 / 160,
48 PI2 / 176,
49 PI2 / 208,
50 PI2 / 240,
51 PI2 / 256,
52 PI2 / 288};
53
54std::string
55TConformalFinder::version(void) const {
56 return "3.03";
57}
58
59TConformalFinder::TConformalFinder(unsigned fastFinder,
60 unsigned slowFinder,
61 unsigned perfectSegmentFinding,
62 float maxSigma,
63 float maxSigmaStereo,
64 float salvageLevel,
65 unsigned minNLinksForSegment,
66 unsigned minNCoreLinks,
67 unsigned minNSegments,
68 unsigned salvageLoadWidth,
69 unsigned stereoMode,
70 unsigned stereoLoadWidth,
71 float szSegmentDistance,
72 float szLinkDistance,
73 unsigned fittingFlag)
74: _fastFinder(fastFinder),
75 _slowFinder(slowFinder),
76 _doT0Reset(false),
77 _perfectSegmentFinding(perfectSegmentFinding),
78 _segmentSeparation(4),
79 _minNLinksForSegment(minNLinksForSegment),
80 _minNLinksForSegmentInRefine(2), //3, 20060307
81 _maxNLinksForSegment(8),
82 _maxWidthForSegment(4),
83 _minNCoreLinks(minNCoreLinks),
84 _minNSegments(minNSegments),
85 _minUsedFractionSlow2D(0.5),
86 _salvageLoadWidth(salvageLoadWidth),
87 _stereoMode(stereoMode),
88 _stereoLoadWidth(stereoLoadWidth),
89 _appendLoad(2),
90 _maxSigma2(maxSigma),
91 _T0ResetDone(false),
92 _builder("conformal builder",
93 maxSigma,
94 maxSigmaStereo,
95 salvageLevel,
96 szSegmentDistance,
97 szLinkDistance,
98 fittingFlag),
99 _s(0)
100#ifdef TRKRECO_WINDOW
101 , _rphiWindow("rphi window")
102#endif
103{
104 _linkMaxDistance[0] = 0.02;
105 _linkMaxDistance[1] = 0.025;
106 _linkMaxDistance[2] = 0.025;
107 _linkMinDirAngle[0] = 0.97; //0.98
108 _linkMinDirAngle[1] = 0.95; //0.97
109 _linkMinDirAngle[2] = 0.95; //0.97
110 // BASF_Sharedmem->allocate("TrkConfSum", sizeof(struct summary));
111
112 //yzhang add 2012-05-03
113 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
114 if(scmgn!=StatusCode::SUCCESS) {
115 std::cout<< __FILE__<<"Unable to open Magnetic field service"<<std::endl;
116 }
117}
118
121
122void
123TConformalFinder::dump(const std::string & msg, const std::string & pre) const {
124 if (msg.find("state") != std::string::npos) {
125 std::cout << pre;
126 TFinderBase::dump(msg);
127 std::cout << pre;
128 std::cout << "#axial=" << _hits[0].length();
129 std::cout << ",#stereo=" << _hits[1].length();
130 }
131 if (msg.find("summary") != std::string::npos || msg.find("detail") != std::string::npos) {
132/* struct summary s;
133 // bzero((char *) & s, sizeof(struct summary));
134 memset((char *) & s, 0, sizeof(struct summary));
135 for (int i = 0; i < BASF_Sharedmem->nprocess(); i++) {
136 int size;
137 struct summary & r = * (struct summary *)
138 BASF_Sharedmem->get_pointer(i, "TrkConfSum", & size);
139 s._nEvents += r._nEvents;
140 s._nTracksFast3D += r._nTracksFast3D;
141 s._nTracksSlow3D += r._nTracksSlow3D;
142 s._nTracksFast2D += r._nTracksFast2D;
143 s._nTracksSlow2D += r._nTracksSlow2D;
144 s._nTracksFast2DBadQuality += r._nTracksFast2DBadQuality;
145 s._nTracksSlow2DBadQuality += r._nTracksSlow2DBadQuality;
146 }
147
148 std::cout << pre
149 << "all events : " << s._nEvents << std::endl;
150 std::cout << pre
151 << "fast 3D tracks : "
152 << s._nTracksFast3D << std::endl;
153 std::cout << pre
154 << "fast 2D tracks(line failed) : "
155 << s._nTracksFast2D << std::endl;
156 std::cout << pre
157 << "fast 2D tracks(3D failed) : "
158 << s._nTracksFast2DBadQuality << std::endl;
159 std::cout << pre
160 << "slow 3D tracks : "
161 << s._nTracksSlow3D << std::endl;
162 std::cout << pre
163 << "slow 2D tracks(line failed) : "
164 << s._nTracksSlow2D << std::endl;
165 std::cout << pre
166 << "slow 2D tracks(3D failed) : "
167 << s._nTracksSlow2DBadQuality << std::endl;
168*/
169 }
170}
171
172void
173TConformalFinder::clear(void) {
174 for (unsigned i = 0; i < 3; i++) {
175 if (i == 2)
176 HepAListDeleteAll(_allHits[i]);
177 else
178 _allHits[i].removeAll();
179 _hits[i].removeAll();
180 _unused[i].removeAll();
181 }
182 for (unsigned i = 0; i < 2; i++) {
183 for (unsigned j = 0; j < 6; j++) {
184#ifdef TRKRECO_DEBUG_DETAIL
185 cout<<"_allSegments length = "<<_allSegments[i][j].length()<<" _allUnused length = "<<_allUnused[i][j].length()<<endl;
186#endif
187 HepAListDeleteAll(_allSegments[i][j]);
188 _allUnused[i][j].removeAll();
189 }
190 }
191 _2DTracks.removeAll();
192 _3DTracks.removeAll();
193}
194
195void
196TConformalFinder::selectGoodHits(void) {
197 for (unsigned i = 0; i < 2; i++) {
198 unsigned n = _allHits[i].length();
199 for (unsigned j = 0; j < n; j++) {
200 unsigned state = _allHits[i][j]->hit()->state();
201// if (_allHits[i][j]->hit()->reccdc()->tdc>800) continue;
202// if ((state & WireHitIsolated) && (state & WireHitContinuous)) //tmp...
203// if ((state & WireHitIsolated))
204 _hits[i].append(_allHits[i][j]);
205// else
206 _unused[i].append(_allHits[i][j]);
207 }
208 }
209 _hits[2].append(_hits[0]);
210 _hits[2].append(_hits[1]);
211 _unused[2].append(_unused[0]);
212 _unused[2].append(_unused[1]);
213
214#ifdef TRKRECO_WINDOW
215 _rphiWindow.clear();
216 _rphiWindow.skip(false);
217 _rphiWindow.skipAllWindow(false);
218 _rphiWindow.append(_allHits[2]);
219 _rphiWindow.append(_hits[2], leda_pink);
220// _rphiWindow.skip(true);
221// _rphiWindow.skipAllWindow(true);
222#endif
223}
224
225void
226TConformalFinder::findSegments(void) {
227
228 //...Create lists of links for each super layer...
229 AList<TMLink> links[11];
230 unsigned n = _hits[2].length();
231 for (unsigned i = 0; i < n; i++) {
232 TMLink & l = * _hits[2][i];
233 links[l.wire()->superLayerId()].append(l);
234 }
235
236 //...Create phi hists...
237 THistogram * hist[11];
238/* hist[0] = new THistogram(64);
239 hist[1] = new THistogram(80);
240 hist[2] = new THistogram(96);
241 hist[3] = new THistogram(128);
242 hist[4] = new THistogram(144);
243 hist[5] = new THistogram(160);
244 hist[6] = new THistogram(176);
245 hist[7] = new THistogram(208);
246 hist[8] = new THistogram(240);
247 hist[9] = new THistogram(256);
248 hist[10] = new THistogram(288);
249*/
250 hist[0] = new THistogram(56);
251 hist[1] = new THistogram(80);
252 hist[2] = new THistogram(88);
253 hist[3] = new THistogram(112);
254 hist[4] = new THistogram(140);
255 hist[5] = new THistogram(160);
256 hist[6] = new THistogram(176);
257 hist[7] = new THistogram(208);
258 hist[8] = new THistogram(240);
259 hist[9] = new THistogram(256);
260 hist[10] = new THistogram(288);
261
262 for (unsigned i = 0; i < 11; i++) {
263// unsigned superLayerId = i / 2;
264// unsigned id = i % 2;
265 unsigned superLayerId;
266 unsigned id;
267 switch(i){
268 case 0: superLayerId = 0; id = 1; break;
269 case 1: superLayerId = 1; id = 1; break;
270 case 2: superLayerId = 0; id = 0; break;
271 case 3: superLayerId = 1; id = 0; break;
272 case 4: superLayerId = 2; id = 0; break;
273 case 5: superLayerId = 2; id = 1; break;
274 case 6: superLayerId = 3; id = 1; break;
275 case 7: superLayerId = 4; id = 1; break;
276 case 8: superLayerId = 5; id = 1; break;
277 case 9: superLayerId = 3; id = 0; break;
278 case 10: superLayerId = 4; id = 0; break;
279 default: break;
280 }
281 // hist[i]->fillX(links[i]);
282 hist[i]->fillPhi(links[i]);
283
284 //...Segment finding...
285 AList<TSegment> tmp = hist[i]->segments();
286
287 //...Remove bad segments...
288 unsigned n = tmp.length();
289 if (id == 0) {
290 AList<TSegment> bad;
291 for (unsigned j = 0; j < n; j++) {
292 TSegment * s = tmp[j];
293 if (s->links().length() < _minNLinksForSegment)
294 bad.append(s);
295 else if (s->links().length() > _maxNLinksForSegment)
296 bad.append(s);
297// else if (Width(s->links()) > _maxWidthForSegment)
298// bad.append(s);
299 }
300 tmp.remove(bad);
301 for (unsigned j = 0; j < bad.length(); j++) {
302 _unused[id].append(bad[j]->links());
303 _unused[2].append(bad[j]->links());
304 }
305 HepAListDeleteAll(bad);
306 n = tmp.length();
307 }
308
309 //...Classify segments...
310// if (n > 1) {
311// for (unsigned k = 0; k < (n - 1); k++) {
312// TSegment & s0 = * tmp[k];
313// bool ok = true;
314// for (unsigned l = k + 1; l < n; l++) {
315// TSegment & s1 = * tmp[l];
316// int distance =
317// abs(InnerMost(s0.links())->wire()->localIdDifference(
318// * InnerMost(s1.links())->wire()));
319// if (distance < _segmentSeparation) {
320// s0.state(s0.state() | TSegmentCrowd);
321// s1.state(s1.state() | TSegmentCrowd);
322// }
323// }
324// }
325// }
326
327 //...Store them...
328 _allSegments[id][superLayerId].append(tmp);
329 _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
330 delete hist[i];
331 }
332
333#ifdef TRKRECO_DEBUG_DETAIL
334 std::cout << "... segment finding finished" << std::endl;
335#endif
336}
337
338void
339TConformalFinder::linkSegments(unsigned level) {
340 //...Clear old links...
341 unsigned superLayer = 5;
342 while (--superLayer) {
343 AList<TSegment> & segments = _allUnused[0][superLayer];
344 unsigned n = segments.length();
345 //cout<<"SL: "<<superLayer<<"; segments.length:"<<n<<endl;
346 for (unsigned i = 0; i < n; i++) {
347 TSegment & base = * segments[i];
348#ifdef TRKRECO_DEBUG
349 if (base.tracks().length()) {
350 std::cout << "TConformalFinder::linkSegments !!! segment has ";
351 std::cout << "an assignment to track(s)" << std::endl;
352 }
353#endif
354 base.innerLinks().removeAll();
355 base.outerLinks().removeAll();
356 }
357 }
358
359 //...Super layer loop...
360 superLayer = 5;
361 while (--superLayer) {
362 AList<TSegment> & segments = _allUnused[0][superLayer];
363 AList<TSegment> & targets = _allUnused[0][superLayer - 1];
364 AList<TSegment> & targets2 = _allUnused[0][superLayer - 2];
365 unsigned n = segments.length();
366 for (unsigned i = 0; i < n; i++) {
367 TSegment & base = * segments[i];
368#ifdef TRKRECO_DEBUG
369 if (base.tracks().length()) {
370 std::cout << "TConformalFinder::linkSegments !!! segment has ";
371 std::cout << "an assignment to track(s)" << std::endl;
372 }
373#endif
374
375 const HepPoint3D & p = base.position();
376 HepVector3D v = base.direction();
377// std::cout<<"superlayer "<<superLayer<<" link length "<<base.outerLinks().length()<<std::endl;
378 if (base.outerLinks().length() == 1)
379 v = p - OuterMostUniqueLink(base)->position();
380// if (base.outerLinks().length() == 1)
381// v = p - base.outerLinks()[0]->position();
382// if (base.outerLinks().length() == 1)
383// v = direction(base);
384
385 AList<TSegment> & candidates = base.innerLinks();
386 TSegment * best = link(base, p, v, targets, candidates, level);
387 if ((best == NULL) && (level > 0) && (superLayer > 1)) {
388 best = link(base, p, v, targets2, candidates, level);
389 }
390 if (best) candidates.insert(best);
391// base.dump("hits");
392// if(best)best->dump("hits");
393
394 unsigned m = candidates.length();
395 for (unsigned j = 0; j < m; j++)
396 candidates[j]->outerLinks().append(base);
397 }
398 }
399
400 //check innerLinks. Liuqg
401/* superLayer = 5;
402 while(--superLayer){
403 //cout<<"SuperLayer: "<<superLayer<<endl;
404 AList<TSegment> & segments = _allUnused[0][superLayer];
405 unsigned n = segments.length();
406 for (unsigned k = 0; k < n; k++) {
407 TSegment & base = * segments[k];
408 int inners = base.innerLinks().length();
409 if (inners<2) continue;
410 AList<TSegment> badList;
411 for(int i = 0; i < inners; ++i){
412 for(int j = i+1; j < inners; ++j){
413 AList<TMLink> links = base.innerLinks()[j]->links();
414 AList<TMLink> multiLinks;
415 for (int k = 0; k < links.length(); ++k){
416 TMLink & l = * links[k];
417 if (base.innerLinks()[i]->links().hasMember(l)) multiLinks.append(l);
418 }
419 if (multiLinks.length() < 2) continue;
420 int nHits[4];
421 int multiLayers = 0;
422 for (int k = 0; k < multiLinks.length(); ++k)
423 ++nHits[multiLinks[k]->hit()->wire()->localLayerId()];
424 for (int k = 0; k < 4; ++k)
425 if (nHits[k] > 0) ++multiLayers;
426 if(multiLayers >= 2 && (links.length() > base.innerLinks()[i]->links().length()))
427 badList.append(base.innerLinks()[i]);
428 else if (multiLayers >= 2) badList.append( base.innerLinks()[j]);
429 else continue;
430 }
431 }
432 base.innerLinks().remove(badList);
433 }
434 }*/
435
436#ifdef TRKRECO_DEBUG
437 superLayer = 5;
438 while(--superLayer){
439 cout<<"SuperLayer: "<<superLayer<<endl;
440 AList<TSegment> & segments = _allUnused[0][superLayer];
441 unsigned n = segments.length();
442 for (unsigned i = 0; i < n; i++) {
443 TSegment & base = * segments[i];
444 cout<<"Layer: "<<base.links()[0]->hit()->wire()->layerId()
445 <<" Local: "<<base.links()[0]->hit()->wire()->localId()
446 <<" innerSeg: "<<base.innerLinks().length()<<endl;
447 }
448 }
449#endif
450
451#ifdef TRKRECO_WINDOW
452 // _rphiWindow.skipAllWindow(false);
453 displayStatus("Conf::linkSegments ... results");
454 _rphiWindow.wait();
455#endif
456}
457
458void
459TConformalFinder::resolveSegments(AList<TTrack> & trackCandidates) const {
460#ifdef TRKRECO_DEBUG
461 std::cout << "... resolving assignments" << std::endl;
462#endif
463
464 //...Pick up segments which has multiple assignments...
465 AList<TSegment> multi;
466 unsigned n = trackCandidates.length();
467 for (unsigned i = 0; i < n; i++) {
468 TTrack & t = * trackCandidates[i];
469 AList<TSegment> & segments = t.segments();
470 unsigned nS = segments.length();
471 for (unsigned j = 0; j < nS; j++) {
472 if (segments[j]->tracks().length() > 1)
473 multi.append(segments[j]);
474 }
475 }
476 multi.purge();
477
478 //...Resolve assignments...
479 n = multi.length();
480 for (unsigned i = 0; i < n; i++) {
481 TSegment & s = * multi[i];
482 const AList<TTrack> & tracks = s.tracks();
483 unsigned nT = tracks.length();
484
485 //...Check TMLink overlap...
486 AList<TMLink> multiLinks;
487 const AList<TMLink> & links = s.links();
488 unsigned nL = links.length();
489 for (unsigned j = 0; j < nL; j++) {
490 TMLink & l = * links[i];
491 bool overlap = false;
492 for (unsigned k = 0; k < nT; k++) {
493 TTrack & t = * tracks[k];
494 if (t.links().hasMember(l))
495 overlap = true;
496 }
497 multiLinks.append(l);
498 }
499 multiLinks.purge();
500
501#ifdef TRKRECO_DEBUG
502 std::cout << " segment : ";
503 s.dump("hits sort flag");
504 std::cout << " # of assigned tracks = " << nT << std::endl;
505 std::cout << " # of overlap TMLinks = " << multiLinks.length();
506 std::cout << std::endl;
507#endif
508 //...Select the closest...
509 nL = multiLinks.length();
510 for (unsigned j = 0; j < nL; j++) {
511 TMLink & l = * multiLinks[j];
512
513 float bestDiff = 999999999.;
514 TTrack * best = NULL;
515 for (unsigned k = 0; k < nT; k++) {
516 TTrack & t = * tracks[k];
517 t.approach(l);
518 float distance = (l.positionOnWire() - l.positionOnTrack())
519 .mag();
520 float diff = fabs(distance - l.hit()->drift());
521 if (diff < bestDiff) {
522 bestDiff = diff;
523 best = & t;
524 }
525 }
526
527 for (unsigned k = 0; k < nT; k++) {
528 TTrack & t = * tracks[k];
529 if (& t == best) continue;
530 t.remove(l);
531 }
532
533#ifdef TRKRECO_DEBUG
534 {
535 std::cout << " " << l.wire()->name() << " -> ";
536 std::cout << best->name() << std::endl;
537 }
538#endif
539 }
540 }
541}
542
543AList<TSegment>
544TConformalFinder::removeBadSegments(TTrack & t) const {
545#ifdef TRKRECO_DEBUG
546 std::cout << "... removing bad segments : # used seg = ";
547 std::cout << t.segments().length() << std::endl;
548#endif
549#ifdef TRKRECO_DEBUG_DETAIL
550 for (unsigned i = 0; i < t.segments().length(); i++)
551 t.segments()[i]->dump("hits sort flag", " ");
552#endif
553
554 const AList<TSegment> & segments = t.segments();
555 AList<TSegment> bads;
556 unsigned used = 0;
557 TSegment * innerMost;
558 unsigned n = segments.length();
559 for (unsigned i = 0; i < n; i++) {
560 TSegment & s = * segments[i];
561 AList<TMLink> links = Links(s, t);
562 unsigned nLinks = links.length();
563
564 unsigned nCores = Cores(links).length();
565 unsigned nCoresSegment = s.nCores();
566// if (nCores < _minNCoreLinks) {
567// if (nCores < (nCoresSegment / 2)) {
568 if ((nCores < _minNCoreLinks) && (nCores < ((nCoresSegment + 1) / 2))){
569 bads.append(s);
570 continue;
571 }
572
573 unsigned id = segments[i]->superLayerId();
574 used |= (1 << id);
575 if (! id) innerMost = & s;
576 }
577
578// //...Check used super layers...
579// n = segments.length();
580// for (unsigned i = 0; i < n; i++)
581// if ((used & 0x155) == 0x101)
582// bads.append(innerMost);
583
584 if (! bads.length()) return bads;
585
586#ifdef TRKRECO_DEBUG
587 std::cout << " bad segments : # = " << bads.length() << std::endl;
588#endif
589
590 n = bads.length();
591 for (unsigned i = 0; i < n; i++) {
592 TSegment & s = * bads[i];
593
594#ifdef TRKRECO_DEBUG
595 AList<TMLink> links = Links(s, t);
596 unsigned nLinks = links.length();
597 unsigned nCores = Cores(links).length();
598 std::cout << " # used links = " << nLinks;
599 std::cout << ", # used cores = " << nCores << std::endl;
600 s.dump("hits sort flag", " ");
601#endif
602
603 s.tracks().remove(t);
604 t.segments().remove(s);
605 t.remove(s.links());
606 }
607
608 //...Refit...
609 t.fit();
610
611#ifdef TRKRECO_DEBUG
612 std::cout << "... refitting" << std::endl;
613 t.dump("detail", "2nd> ");
614#endif
615
616 return bads;
617}
618
619TSegment *
620TConformalFinder::link(const TSegment & base,
621 const HepPoint3D & p,
622 const HepVector3D & v,
623 const AList<TSegment> & candidates,
624 AList<TSegment> & alternatives,
625 unsigned level) const {
626
627#ifdef TRKRECO_DEBUG
628 std::cout << " link : base = " << std::endl;
629 base.dump("vector hits mc sort", "-> ");
630 std::cout << " p=" << p << ", v=" << v.unit() << std::endl;
631#endif
632
633 //...Parameters...
634 //static const float maxDistance = 0.02;
635 //static const float minDirAngle = 0.97;
636
637 //...Candidate loop...
638 unsigned n = candidates.length();
639 float maxAngle = -999.;
640 float minDistance = _linkMaxDistance[level];
641 float maxDirAngle = _linkMinDirAngle[level];
642
643 TSegment * best = NULL;
644 for (unsigned j = 0; j < n; j++) {
645 TSegment & current = * candidates[j];
646
647#ifdef TRKRECO_DEBUG_DETAIL
648 current.dump("vector hits mc sort", " ");
649 std::cout << " dist,dirAngle,angle=" << current.distance(p, v);
650 std::cout << "," << ((current.position() - p).unit()).dot(v.unit());
651 std::cout << "," << v.dot(current.direction()) << std::endl;
652#endif
653 float distance = current.distance(p, v);
654 //cout<<"current, Layer: "<<current.links()[0]->hit()->wire()->layerId()
655 // <<" local: "<<current.links()[0]->hit()->wire()->localId()<<endl;
656 //cout<<" distance of segs = "<<distance<<" limit = "<<_linkMaxDistance[level]<<endl;
657 if (distance > _linkMaxDistance[level]) continue;
658
659 HepVector3D dir = (current.position() - p);
660 if (dir.x() > M_PI) dir.setX(dir.x() - PI2);
661 else if (dir.x() < - M_PI) dir.setX(PI2 + dir.x());
662
663 float dirAngle = dir.unit().dot(v.unit());
664 //cout<<" dirAngle of segs = "<<dirAngle<<" limit = "<<_linkMinDirAngle[level]<<endl;
665 if (dirAngle < _linkMinDirAngle[level]) continue; //tmp for tsf
666
667// if (distance < minDistance) {
668// minDistance = distance;
669// best = & current;
670// }
671 if (dirAngle > maxDirAngle) {
672 maxDirAngle = dirAngle;
673 best = & current;
674 }
675 alternatives.append(current);
676 }
677
678#ifdef TRKRECO_DEBUG
679 std::cout << " # of best + alternatives = " << alternatives.length() << std::endl;
680 if (best) {
681 std::cout << " Best is ";
682 best->dump("hits");
683 }
684#endif
685
686 alternatives.remove(best);
687 if (best) return best;
688 return NULL;
689}
690
691void
692TConformalFinder::deleteTrack(TTrack & t) const {
693 const AList<TSegment> & segments = t.segments();
694 unsigned n = segments.length();
695 for (unsigned i = 0; i < n; i++) {
696 TSegment & s = * segments[i];
697 s.tracks().remove(t);
698 }
699
700 delete & t;
701}
702
703void
704TConformalFinder::removeUsedSegments(const AList<TTrack> & tracks) {
705 unsigned n = tracks.length();
706 for (unsigned i = 0; i < n; i++) {
707 TTrack & t = * tracks[i];
708 AList<TSegment> & segments = t.segments();
709 AList<TSegment> toBeRemoved;
710 unsigned nS = segments.length();
711 for (unsigned j = 0; j < nS; j++) {
712 TSegment & s = * segments[j];
713 unsigned sId = s.superLayerId();
714// unsigned as = sId % 2;
715// unsigned id = sId / 2;
716 unsigned as;
717 unsigned id;
718 switch (sId) {
719 case 0: as = 1; id = 0;break;
720 case 1: as = 1; id = 1;break;
721 case 2: as = 0; id = 0;break;
722 case 3: as = 0; id = 1;break;
723 case 4: as = 0; id = 2;break;
724 case 5: as = 1; id = 2;break;
725 case 6: as = 1; id = 3;break;
726 case 7: as = 1; id = 4;break;
727 case 8: as = 1; id = 5;break;
728 case 9: as = 0; id = 3;break;
729 case 10: as = 0; id = 4;break;
730 default: break;
731 }
732
733 //...Check used links...
734 AList<TMLink> links = Links(s, t);
735 if (links.length() == 0) {
736 s.tracks().remove(t);
737 toBeRemoved.append(s);
738#ifdef TRKRECO_DEBUG
739 std::cout << "!!! why this happends???" << std::endl;
740 std::cout << " no used link" << std::endl;
741#endif
742 continue;
743 }
744
745 //...Remove from lists...
746 _allUnused[as][id].remove(s);
747
748 //...Remove incoming links...
749 const AList<TSegment> & outers = s.outerLinks();
750 unsigned nO = outers.length();
751// std::cout << " >>> removing outer links" << std::endl;
752// s.dump("hits", " To ");
753 for (unsigned i = 0; i < nO; i++) {
754 outers[i]->innerLinks().remove(s);
755// outers[i]->dump("hits", " From ");
756 }
757
758 //...Remaining hits...
759 AList<TMLink> unused = s.links();
760 unused.remove(links);
761 s.remove(unused);
762 unsigned nUnused = unused.length();
763
764 //...Create new segment if too many remaining links...
765/* if (nUnused < _minNLinksForSegment) {
766 for (unsigned k = 0; k < nUnused; k++) {
767 _unused[as].append(unused[k]);
768 _unused[2].append(unused[k]);
769 }
770 }
771 else {
772 TSegment * ss = new TSegment(unused);
773 AList<TSegment> newSegments = ss->split(); //tmp for tsf
774 if (newSegments.length() == 0) {
775 ss->solveDualHits();
776 newSegments.append(ss);
777 for (unsigned k = 0; k < nO; k++)
778 outers[k]->innerLinks().append(ss);
779 }
780 else {
781 delete ss;
782 }
783 _allUnused[as][id].append(newSegments);
784 _allSegments[as][id].append(newSegments);
785 }
786*/
787 }
788 segments.remove(toBeRemoved);
789 }
790}
791
792void
793TConformalFinder::updateTLinks(AList<TTrack> & tracks) {
794 unsigned n = tracks.length();
795 for (unsigned i = 0; i < n; i++) {
796 const AList<TMLink> & links = tracks[i]->links();
797 unsigned nL = links.length();
798 for (unsigned j = 0; j < nL; j++) {
799 TMLink & l = * links[j];
800 tracks[i]->approach(l);
801 }
802 }
803}
804
805int
806TConformalFinder::doit(const AList<TMDCWireHit> & axial,
807 const AList<TMDCWireHit> & stereo,
808 AList<TTrack> & tracks,
809 AList<TTrack> & tracks2D) {
810#ifdef TRKRECO_DEBUG
811 _stage = ConformalInitialization;
812#endif
813
814 static bool first = true;
815 if (first) {
816 first = false;
817 int size;
818// _s = (struct summary *)
819// BASF_Sharedmem->get_pointer(BASF_Sharedmem->get_id(),
820// "TrkConfSum",
821// & size);
822 }
823
824 //...For debug...
825 if (debugLevel()) {
826 std::cout << name() << " ... processing"
827 << " axial=" << axial.length()
828 << ",stereo=" << stereo.length()
829 << ",tracks=" << tracks.length()
830 << std::endl;
831 }
832
833 //...Create TMLinks with conformal position...
834// TConformalFinder0::conformalTransformation(ORIGIN, axial, _allHits[0]);
835// TConformalFinder0::conformalTransformation(ORIGIN, stereo, _allHits[1]);
836 // for Tsf
837 TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, axial, _allHits[0]);
838 TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, stereo, _allHits[1]);
839 _allHits[2].append(_allHits[0]);
840 _allHits[2].append(_allHits[1]);
841
842 //...Select good hits...
843 selectGoodHits();
844#ifdef TRKRECO_DEBUG
845 cout<<"axial Hits:"<<_allHits[0].length()<<" good axial hits:"<<_hits[0].length()
846 <<endl<<"stereo Hits:"<<_allHits[1].length()<<" good setero hits:"<<_hits[1].length()
847 <<endl;
848#endif
849
850 //...Segment finding...
851 if (_perfectSegmentFinding)
852 findSegmentsPerfect();
853 else
854// findSegments();
855 findSegmentsTsf();
856
857#ifdef TRKRECO_DEBUG
858 cout<<"axial Seg: "<<_allSegments[0][0].length()<<", "<<_allSegments[0][1].length()
859 <<", "<<_allSegments[0][2].length()<<", "<<_allSegments[0][3].length()
860 <<", "<<_allSegments[0][4].length()<<", "<<_allSegments[0][5].length()
861 <<" stereo Seg: "<<_allSegments[1][0].length()<<", "<<_allSegments[1][1].length()
862 <<", "<<_allSegments[1][2].length()<<", "<<_allSegments[1][3].length()
863 <<", "<<_allSegments[1][4].length()<<", "<<_allSegments[1][5].length()<<endl;
864
865 for(int i=0; i<2; ++i) {
866 if (i == 0) cout<<"Axial......"<<endl;
867 else cout<<"Stereo......"<<endl;
868 for(int j=0; j<6; ++j) {
869 for(int k=0; k<_allSegments[i][j].length(); ++k){
870 cout<<"SL(a/s): "<<j<<" seeds in Seg"<<k<<": "<<_allSegments[i][j][k]->links().length()<<endl;
871 for(int kk=0; kk<_allSegments[i][j][k]->links().length(); ++kk)
872 cout<<" layerId: "<<_allSegments[i][j][k]->links()[kk]->hit()->wire()->layerId()
873 <<" localId: "<<_allSegments[i][j][k]->links()[kk]->hit()->wire()->localId()<<endl;
874 }
875 }
876 }
877#endif
878
879 //...Fast finding...
880 unsigned level = 0;
881 _T0ResetDone = false;
882 if (_fastFinder) {
883 linkSegments(level);
884 fastFinding2D(level);
885 updateTLinks(_2DTracks);
886
887 //...T0 reset here...
888 if (_doT0Reset) {
889 std::cout << "TConformalFinder ... T0 reset is done" << std::endl;
890 _T0ResetDone = true;
891 return 0;
892 }
893
894#ifdef TRKRECO_WINDOW
895 _rphiWindow.skip(false);
896#endif
897 fastFinding3D(level);
898 updateTLinks(_2DTracks);
899 updateTLinks(_3DTracks);
900
901#ifdef TRKRECO_WINDOW
902 _rphiWindow.skip(false);
903 displayTracks(_2DTracks, _allUnused, "all 2D after fast level 0");
904 displayTracks(_3DTracks, _allUnused, "all 3D after fast level 0");
905#endif // TRKRECO_WINDOW
906
907
908 //...Fast finding again...
909 level = 1;
910 linkSegments(level);
911 fastFinding2D(level);
912 updateTLinks(_2DTracks);
913
914#ifdef TRKRECO_WINDOW
915 _rphiWindow.skip(false);
916#endif
917 fastFinding3D(level);
918 updateTLinks(_2DTracks);
919 updateTLinks(_3DTracks);
920
921#ifdef TRKRECO_WINDOW
922 _rphiWindow.skip(false);
923 displayTracks(_2DTracks, _allUnused, "all 2D after fast level 1");
924 displayTracks(_3DTracks, _allUnused, "all 3D after fast level 1");
925#endif //TRKRECO_WINDOW
926 }
927
928 //...Slow finding...
929 if (_slowFinder) {
930 level = 2;
931 linkSegments(level);
932 slowFinding2D(level);
933 updateTLinks(_2DTracks);
934#ifdef TRKRECO_WINDOW
935 _rphiWindow.skip(false);
936// _rphiWindow.skipAllWindow(false);
937#endif
938 fastFinding3D(level);
939 updateTLinks(_2DTracks);
940 updateTLinks(_3DTracks);
941
942#ifdef TRKRECO_WINDOW
943 _rphiWindow.skip(false);
944// _rphiWindow.skipAllWindow(false);
945 displayTracks(_2DTracks, _allUnused, "all 2D after slow level 2");
946 displayTracks(_3DTracks, _allUnused, "all 3D after slow level 2");
947#endif //TRKRECO_WINDOW
948 }
949
950#ifdef TRKRECO_DEBUG
951 int zsltrk = _2DTracks.length();
952 std::cout << name() << " ... # 3D tracks = " << _3DTracks.length()
953 << ", # 2D tracks = " << _2DTracks.length() << std::endl;
954 for (unsigned i = 0; i < zsltrk; i++){
955// _2DTracks[i]->dump();
956 cout<<"pt:"<<_2DTracks[i]->pt()<<" impact:"<<_2DTracks[i]->impact()<<endl;
957 // nTrack->fill(0., 1.0);
958 }
959#endif
960
961 //...Mask hits with bad chisq...
962 TTrackManager::maskBadHits(_3DTracks, _maxSigma2);
963
964 tracks2D = _2DTracks;
965 tracks = _3DTracks;
966 for(int i=0;i<_3DTracks.length();i++){
967 _3DTracks[i]->finder(TrackOldConformalFinder);
968 _3DTracks[i]->setFinderType(2);
969 }
970#ifdef TRKRECO_DEBUG
971// if (_3DTracks.length() == 0) {
972// cout<<"3D failed: 2D length = "<<_2DTracks.length()<<endl;
973 for(int nn = 0; nn < tracks2D.length(); ++nn){
974 cout<<"2D: "<<nn<<" radius: "<<tracks2D[nn]->radius()<<endl;
975 for (int mm = 0; mm < tracks2D[nn]->links().length(); ++mm)
976 cout<<"layer: "<<tracks2D[nn]->links()[mm]->wire()->layerId()
977 <<" local: "<<tracks2D[nn]->links()[mm]->wire()->localId()<<endl;
978 }
979
980 cout<<"unused axial Seg: "<<_allUnused[0][0].length()<<", "<<_allUnused[0][1].length()
981 <<", "<<_allUnused[0][2].length()<<", "<<_allUnused[0][3].length()
982 <<", "<<_allUnused[0][4].length()<<", "<<_allUnused[0][5].length()
983 <<" unused stereo Seg: "<<_allUnused[1][0].length()<<", "<<_allUnused[1][1].length()
984 <<", "<<_allUnused[1][2].length()<<", "<<_allUnused[1][3].length()
985 <<", "<<_allUnused[1][4].length()<<", "<<_allUnused[1][5].length()
986 <<endl;
987
988 for(int i=0; i<2; ++i) {
989 if (i == 0) cout<<"unused axial hits: "<<endl;
990 else cout<<"unused stereo hits: "<<endl;
991 for(int k = 0; k < _allHits[i].length(); ++k) {
992 if (_allHits[i][k]->hit()->state() & WireHitUsed) continue;
993 else cout<<" layerId: "<<_allHits[i][k]->hit()->wire()->layerId()
994 <<" localId: "<<_allHits[i][k]->hit()->wire()->localId()<<endl;
995 }
996 if (i == 0) cout<<"unused axial link in segs: "<<endl;
997 else cout<<"unused stereo link in segs: "<<endl;
998 for(int j=0; j<6; ++j) {
999 for(int k=0; k<_allUnused[i][j].length(); ++k){
1000 cout<<"sl: "<<i<<" "<<j<<" length of seg "<<k<<": "<<_allUnused[i][j][k]->links().length()<<endl;
1001 for(int kk=0; kk<_allUnused[i][j][k]->links().length(); ++kk)
1002 cout<<" layerId: "<<_allUnused[i][j][k]->links()[kk]->hit()->wire()->layerId()
1003 <<" localId: "<<_allUnused[i][j][k]->links()[kk]->hit()->wire()->localId()<<endl;
1004 }
1005 }
1006 }
1007// }
1008#endif
1009
1010 //...Termination...
1011/* tracks = _3DTracks;
1012 tracks2D = _2DTracks;
1013 ++_s->_nEvents;
1014 unsigned n3 = _3DTracks.length();
1015 for (unsigned i = 0; i < n3; i++)
1016 if (_3DTracks[i]->finder() & TrackFastFinder)
1017 ++_s->_nTracksFast3D;
1018 else if (_3DTracks[i]->finder() & TrackSlowFinder)
1019 ++_s->_nTracksSlow3D;
1020 unsigned n2 = _2DTracks.length();
1021 for (unsigned i = 0; i < n2; i++)
1022 if (_2DTracks[i]->finder() & TrackFastFinder)
1023 ++_s->_nTracksFast2D;
1024 else if (_2DTracks[i]->finder() & TrackSlowFinder)
1025 ++_s->_nTracksSlow2D;
1026
1027 if (debugLevel()) {
1028 std::cout << name() << " ... # 3D tracks = " << _3DTracks.length()
1029 << ", # 2D tracks = " << _2DTracks.length() << std::endl;
1030 }
1031
1032#ifdef TRKRECO_DEBUG
1033 _stage = ConformalOutside;
1034#endif
1035*/
1036 return 0;
1037}
1038
1039void
1040TConformalFinder::fastFinding3D(unsigned level) {
1041
1042#ifdef TRKRECO_DEBUG_DETAIL
1043 std::cout << "... fast finding (3D) : level = " << level << std::endl;
1044#endif
1045#ifdef TRKRECO_DEBUG
1046 if (_stage == ConformalSlow2D) {
1047 _stage = ConformalSlow3D;
1048 }
1049 else {
1050 _stage = ConformalFast3DLevel0;
1051 if (level == 1)
1052 _stage = ConformalFast3DLevel1;
1053 }
1054#endif
1055
1056 AList<TTrack> tracks3D;
1057 AList<TTrack> touched;
1058 AList<TSegment> bads;
1059 unsigned n = _2DTracks.length();
1060 for (unsigned i = 0; i < n; i++) {
1061
1062 const TTrack & t = * _2DTracks[i];
1063
1064#ifdef TRKRECO_DEBUG
1065 cout<<"links in 2Dtracks: "<<t.links().length()<<endl;
1066 for (int ii = 0; ii < t.links().length(); ++ii) {
1067 cout<<"layerId: "<<t.links()[ii]->wire()->layerId()
1068 <<" localId: "<<t.links()[ii]->wire()->localId()<<endl;
1069 }
1070 cout<<"center: "<<t.center()<<" radius = "<<t.radius()<<endl;
1071#endif
1072
1073#ifdef TRKRECO_DEBUG_DETAIL
1074 std::cout << "==> fast 3D building : " << t.name() << std::endl;
1075 t.dump("hits sort flag", " 2D track hits=");
1076#endif
1077 AList<TSegment> segments = stereoSegments(t);
1078 //jialk
1079 unsigned NS = segments.length();
1080 if ( NS >= 30 ) continue;
1081
1082 //yuany
1083 //cout<<"stereo segments:"<<segments.length()<<endl;
1084 //unsigned n = segments.length();
1085 //for (unsigned i = 0; i < n; i++) {
1086 //segments[i]->dump("hits");
1087 //}
1088
1089 AList<TSegment> badSegments;
1090 if (_stereoMode == 2)
1091 badSegments = stereoSegmentsFromBadHits(t);
1092 //yuany
1093 //cout<<"bad stereo segments:"<<endl;
1094 //n = badSegments.length();
1095 //for (unsigned i = 0; i < n; i++) {
1096 //badSegments[i]->dump("hits");
1097 //}
1098
1099#ifdef TRKRECO_WINDOW
1100 displayStatus("Conf::fast3D ... seed");
1101 _rphiWindow.append(segments, leda_blue);
1102 _rphiWindow.append(badSegments, leda_red);
1103 _rphiWindow.oneShot(t, leda_green);
1104#endif
1105
1106 //cout<<"stereo mode: "<<_stereoMode<<endl;
1107
1108 //...Save a 2D track...
1109 TTrack * s = NULL;
1110 if (_stereoMode)
1111 s = _builder.buildStereoNew(t, segments, badSegments);
1112 else
1113 s = _builder.buildStereo(t, segments);
1114 HepAListDeleteAll(badSegments);
1115 //cout<<"after 3D build..."<<endl;
1116 if (! s) {
1117#ifdef TRKRECO_DEBUG_DETAIL
1118 std::cout << "... 3D failure" << std::endl;
1119#endif
1120#ifdef TRKRECO_WINDOW
1121 displayStatus("Conf::fastd3D ... 3D failed");
1122 _rphiWindow.append(segments, leda_blue);
1123 _rphiWindow.oneShot(t, leda_red);
1124#endif
1125 continue;
1126 }
1127 //...Quality check...
1128 if (! TTrackManager::goodTrack(* s)) {
1129#ifdef TRKRECO_DEBUG_DETAIL
1130 std::cout << "... 3D failure (bad quality)" << std::endl;
1131#endif
1132#ifdef TRKRECO_WINDOW
1133 displayStatus("Conf::fastd3D ... 3D failed (bad quality)");
1134 _rphiWindow.append(segments, leda_blue);
1135 _rphiWindow.oneShot(* s, leda_red);
1136#endif
1137// Remove following sentences because _s is NULL now
1138// and it's wrong to use the pointer for some noise events... 05/10/15 zangsl
1139//
1140// if (s->finder() & TrackFastFinder)
1141// ++_s->_nTracksFast2DBadQuality;
1142// else if (s->finder() & TrackSlowFinder)
1143// ++_s->_nTracksSlow2DBadQuality;
1144
1145 deleteTrack(* s);
1146 continue;
1147 }
1148
1149 //...New name...
1150 s->name(t.name() + "-3D");
1151
1152#ifdef TRKRECO_WINDOW
1153 displayStatus("Conf::fastd3D ... 3D ok");
1154 _rphiWindow.append(segments, leda_blue);
1155 _rphiWindow.oneShot(* s, leda_green);
1156#endif
1157
1158 //...Salvage by segments...
1159 salvage(* s, 3, bads);
1160 //cout<<"3D finder done!"<<endl;
1161 tracks3D.append(s);
1162 touched.append(_2DTracks[i]);
1163
1164 s->assign(WireHitConformalFinder);
1165 s->quality(0);
1166
1167 //...Segment...
1168 static AList<TTrack> tmp;
1169 tmp.removeAll();
1170 tmp.append(s);
1171 removeUsedSegments(tmp);
1172 // liucy :add constraint of number of tracks
1173 if(tracks3D.length()>20)break;
1174#ifdef TRKRECO_DEBUG_DETAIL
1175 std::cout << "... 3D finished : " << s->name() << std::endl;
1176 s->dump("detail", " ");
1177#endif
1178#ifdef TRKRECO_WINDOW
1179 displayStatus("Conf::fastd3D ... finished");
1180 _rphiWindow.oneShot(* s, leda_green);
1181#endif
1182 }
1183
1184 // resolveSegments(tracks3D);
1185 _3DTracks.append(tracks3D);
1186 _2DTracks.remove(tracks3D);
1187 _2DTracks.remove(touched);
1188 for (unsigned i = 0; i < touched.length(); i++) {
1189 deleteTrack(* touched[i]);
1190 // saved[i]->fit();
1191 }
1192}
1193
1194#ifdef TRKRECO_WINDOW
1195void
1196TConformalFinder::displayStatus(const std::string & m) const {
1197 _rphiWindow.clear();
1198 _rphiWindow.text(m);
1199 _rphiWindow.append(_allHits[2], leda_pink);
1200 _rphiWindow.append(_hits[2], leda_cyan);
1201 _rphiWindow.append(_unused[2]);
1202 displayAppendSegments(_allSegments, leda_grey1);
1203 displayAppendSegments(_allUnused, leda_orange);
1204}
1205
1206void
1207TConformalFinder::displayAppendSegments(const AList<TSegment> a[2][6],
1208 leda_color c) const {
1209 for (unsigned i = 0; i < 2; i++) {
1210 for (unsigned j = 0; j < 6; j++) {
1211 const AList<TSegment> & segments = a[i][j];
1212 unsigned n = segments.length();
1213 for (unsigned k = 0; k < n; k++) {
1214 _rphiWindow.append(* segments[k], c);
1215 }
1216 }
1217 }
1218}
1219
1220void
1221TConformalFinder::displayTracks(const AList<TTrack> & l,
1222 const AList<TSegment> segments[2][6],
1223 const std::string & c) const {
1224 displayStatus("Conf::display ... " + c);
1225 for (unsigned i = 0; i < l.length(); i++)
1226 _rphiWindow.append(* l[i], leda_green);
1227 _rphiWindow.wait();
1228}
1229#endif
1230
1231bool
1232TConformalFinder::quality2D(TTrack & t, unsigned level) const {
1233#ifdef TRKRECO_WINDOW
1234 std::string s = "Conf::quality2D ... level " + itostring(level);
1235#endif
1236 //jialk caution origin is 15.
1237 //...Check # radius...Liuqg
1238 if(fabs(t.radius()) < 15.) return false; //corresponding to pt = 300*1T*0.50/2(m).
1239
1240 //...Check # of segments(superlayers)...
1241 // unsigned nSuperLayers = NSuperLayers(t.links(), 2);
1242 unsigned nSuperLayers = NSuperLayers(t.links());
1243 if (nSuperLayers < _minNSegments) {
1244#ifdef TRKRECO_DEBUG
1245 std::cout << "... quality check : level=" << level << " : bad" << std::endl;
1246 std::cout << " reason : # segments(superlayer) =" << nSuperLayers;
1247 std::cout << " < " << _minNSegments << std::endl;
1248#endif
1249#ifdef TRKRECO_WINDOW
1250 s += " rejected because of few segments(superlayers)";
1251 displayStatus(s);
1252 _rphiWindow.oneShot(t, leda_red);
1253#endif
1254 return false;
1255 }
1256 if (level == 0) {
1257#ifdef TRKRECO_WINDOW
1258 s += " ok : # of used segments(superlayers) = ";
1259 s += itostring(nSuperLayers);
1260 s += " > " + itostring(_minNSegments);
1261 displayStatus(s);
1262 _rphiWindow.oneShot(t, leda_green);
1263#endif
1264
1265 return true;
1266 }
1267
1268 //...Check superlayer usage...
1269 unsigned n3 = NSuperLayers(t.links(), 3);
1270
1271 //...Check super layer usage...
1272 if (n3 != nSuperLayers) {
1273 unsigned sl = SuperLayer(t.links(), 2);
1274 // unsigned sl = SuperLayer(t.links(), 3);
1275 if (sl == 0) {
1276#ifdef TRKRECO_DEBUG
1277 std::cout << "... quality check : level = " << level << " : bad";
1278 std::cout << std::endl;
1279 std::cout << " reason : super layer pattern(n2) = 0 " << std::endl;
1280#endif
1281#ifdef TRKRECO_WINDOW
1282 s += " rejected because of bad super-layer pattern(n2=0)";
1283 displayStatus(s);
1284 _rphiWindow.oneShot(t, leda_red);
1285#endif
1286 return false;
1287 }
1288 //unsigned fl = 0; //belle, liuqg
1289 //while ((sl & (1 << (fl * 2))) == 0) ++fl;
1290 //bool empty = false;
1291
1292 /*unsigned fl[5] = {2, 3, 4, 9, 10}; //tmp....liuqg
1293 unsigned nfl = 0;
1294 for (unsigned i = 0; i < 5; i++) {
1295 if ((sl & (1 << fl[i])) == 1) {
1296 nfl = i;
1297 break;
1298 }
1299 }
1300 bool empty = false;
1301 for (unsigned i = nfl + 1; i < 5; i++) {
1302 bool thisLayer = (sl & (1 << fl[i]));
1303 if (thisLayer && empty) {
1304
1305//#ifdef TRKRECO_DEBUG
1306 std::cout << "... quality check : level = " << level;
1307 std::cout << " : bad" << std::endl;
1308 std::cout << " reason : super layer pattern = ";
1309// for (unsigned j = 0; j < 6; j++) std::cout << (sl >> (j * 2)) % 2;
1310 std::cout << std::endl;
1311//#endif
1312#ifdef TRKRECO_WINDOW
1313 s += " rejected because of bad super-layer pattern";
1314 displayStatus(s);
1315 _rphiWindow.oneShot(t, leda_red);
1316#endif
1317 return false;
1318 }
1319 if (! thisLayer) empty = true;
1320 }*/
1321 }
1322
1323#ifdef TRKRECO_DEBUG
1324 std::cout << "... quality check : level = " << level;
1325 std::cout << " : ok, super layer pattern = " << std::endl;
1326 std::string ptn;
1327 unsigned sl = SuperLayer(t.links(), 2);
1328 for (unsigned j = 0; j < 6; j++) {
1329 unsigned k = (sl >> (j * 2)) % 2;
1330 std::cout << k;
1331 if (k) ptn += "1";
1332 else ptn += "0";
1333 }
1334 std::cout << std::endl;
1335#endif
1336#ifdef TRKRECO_WINDOW
1337 s += " ok : super layer ptn = " + ptn;
1338 displayStatus(s);
1339 _rphiWindow.oneShot(t, leda_green);
1340#endif
1341
1342 return true;
1343}
1344
1345void
1346TConformalFinder::fastFinding2D(unsigned level) {
1347
1348#ifdef TRKRECO_DEBUG_DETAIL
1349 std::cout << "... fast finding (2D) : level = " << level << std::endl;
1350#endif
1351#ifdef TRKRECO_DEBUG
1352 _stage = ConformalFast2DLevel0;
1353 if (level == 1)
1354 _stage = ConformalFast2DLevel1;
1355#endif
1356
1357 AList<TSegment> * segments = & _allUnused[0][0];
1358
1359 unsigned idBase = _2DTracks.length() + _3DTracks.length();
1360 unsigned nTrial = 0;
1361 unsigned seedLayer = 5;
1362 //yuany
1363/* cout<<"before loop "<<endl;
1364 for (unsigned i = 0; i < 5; i++) {
1365 for(unsigned j=0;j<_allUnused[0][i].length();j++){
1366 (_allUnused[0][i])[j]->dump("hits");
1367 (_allUnused[0][i])[j]->dump("link");
1368 }
1369 }*/
1370
1371// while ((seedLayer--) > 2) { 2006/04/27
1372 while ((seedLayer--) > 1) {
1373
1374 //...Seed loop...
1375 AList<TTrack> trackCandidates;
1376 unsigned n = segments[seedLayer].length();
1377
1378 for (unsigned i = 0; i < n; i++) {
1379 TSegment & seed = * segments[seedLayer][i];
1380 std::string name = "f" + itostring(nTrial + idBase);
1381 //yuany
1382 //seed.dump("hits");
1383 //seed.dump("link");
1384#ifdef TRKRECO_DEBUG
1385 std::cout << ">>> fast 2D : super layer " << seedLayer << ", ";
1386 std::cout << i << " : ";
1387 seed.dump("link");
1388#endif
1389
1390 //...Check uniqueness...
1391 AList<TSegment> seeds;
1392 //cout<<"NUiqueLinks:"<<NUniqueLinks(seed)<<"; NMajorLinks:"<<NMajorLinks(seed)<<endl;
1393 //cout<<"innerLinks size:"<<seed.innerLinks().length()<<endl;
1394 if (NUniqueLinks(seed) > 0) {
1395 seeds = UniqueLinks(seed);
1396 }
1397 else if (NMajorLinks(seed) > 0) {
1398 seeds = MajorLinks(seed);
1399 }
1400 else {
1401 continue;
1402 }
1403 seeds.append(seed);
1404
1405#ifdef TRKRECO_DEBUG
1406 cout<<"seeds in fastFinding2D: "<<seeds.length()<<endl;
1407 for (int kk = 0; kk < seeds.length(); ++kk)
1408 cout<<"LayerId: "<<seeds[kk]->links()[0]->wire()->layerId()
1409 <<" LocalId: "<<seeds[kk]->links()[0]->wire()->localId()<<endl;
1410#endif
1411
1412 //...Refine...
1413 refine:
1414
1415 //solve LR of seeds;
1416 for (int k = 0; k < seeds.length(); ++k) //tsf
1417 seeds[k]->solveLR();
1418
1419#ifdef TRKRECO_WINDOW
1420 displayStatus("Conf::fast2D ... seed segments");
1421 _rphiWindow.oneShot(seeds, leda_green);
1422#endif
1423
1424 //...Try to build a track...
1425#ifdef TRKRECO_DEBUG
1426 std::cout << "==> fast building 2D : seed layer = " << seedLayer;
1427 std::cout << ", " << name << std::endl;
1428#endif
1429
1430 TTrack * t = _builder.buildRphi(seeds);
1431 ++nTrial;
1432 //cout<<"buildRphi done."<<endl;
1433
1434 //...Track check...
1435 if (! t) continue;
1436 t->name(name);
1437 bool ok = quality2D(* t, 0);
1438 //cout<<"track quality " <<ok<<endl;
1439 if (! ok) {
1440 deleteTrack(* t);
1441 continue;
1442 }
1443 //cout<<"length of seg in 2d build = "<<t->segments().length()<<" cores in 2d = "<<t->nCores()<<endl;
1444
1445 //...Bad segment rejection...
1446 AList<TSegment> bads = removeBadSegments(* t);
1447 //cout<<"test point no0; bad segment length: "<<bads.length()<<endl;
1448
1449 // This was taken out because the goto loop continues forever for some noise events. 2005/10/14
1450 if (bads.length()) {
1451 ok = quality2D(* t, 1);
1452 if (! ok) {
1453 seeds = refineSegments(* t);
1454 if (seeds.length() >= _minNSegments) {
1455 deleteTrack(* t);
1456 goto refine;
1457 }
1458 deleteTrack(* t);
1459 continue;
1460 }
1461 }
1462 //...Salvage by segments...
1463 salvage(* t, 1, bads);
1464 //cout <<"test point no1"<<endl;
1465 //refineLinks(* t, 3); //change the MinN from 3 to 2, zangshilei 20060307
1466 refineLinks(* t, 2);
1467 //cout<<"test point no1.5"<<endl;
1468
1469 //change the level from 2 to 0, 20060428; because the program stops(MCv5.1.0/run10237.dat)
1470 //release by Liuqg
1471 ok = quality2D(* t, 2);
1472 if (! ok) {
1473 deleteTrack(* t);
1474 continue;
1475 }
1476 //cout<<"test point no2"<<endl;
1477
1478 //...Append segments...
1479 if (NSuperLayers(t->links()) < 5) {
1480 t = expand(* t);
1481 }
1482 //cout<<"test point no3"<<endl;
1483
1484 //t->dump("detail", " ");
1485#ifdef TRKRECO_DEBUG
1486 std::cout << "... 2D finished : " << t->name() << std::endl;
1487 t->dump("hits sort flag pull", " ");
1488#endif
1489#ifdef TRKRECO_WINDOW
1490 displayStatus("Conf::fast2D ... finished");
1491 _rphiWindow.oneShot(* t, leda_green);
1492#endif
1493 t->finder(TrackFastFinder);
1494 t->quality(TrackQuality2D);
1495 t->fitting(TrackFitGlobal);
1496 trackCandidates.append(t);
1497 }
1498
1499 //...Resolve multi-track candidates...
1500 resolveSegments(trackCandidates);
1501
1502 //...Remove used segments...
1503 //yuany
1504/* cout<<"before remove "<<endl;
1505 for (unsigned i = 0; i < 5; i++) {
1506 for(unsigned j=0;j<_allUnused[0][i].length();j++){
1507 (_allUnused[0][i])[j]->dump("hits");
1508 (_allUnused[0][i])[j]->dump("link");
1509 }
1510 }
1511*/
1512 removeUsedSegments(trackCandidates);
1513 //yuany
1514// cout<<"after remove "<<endl;
1515/* for (unsigned i = 0; i < 5; i++) {
1516 for(unsigned j=0;j<_allUnused[0][i].length();j++){
1517 (_allUnused[0][i])[j]->dump("hits");
1518 (_allUnused[0][i])[j]->dump("link");
1519 }
1520 }
1521*/
1522 _2DTracks.append(trackCandidates);
1523 }
1524
1525 resolveHits(_2DTracks);
1526}
1527
1528const TMDCWire *
1529TConformalFinder::conformal2Wire(const HepPoint3D & p) {
1530 static const TMDC & cdc = * TMDC::getTMDC();
1531
1532 std::cout << "p = " << p << std::endl;
1533 float r = sqrt(4. / p.mag2());
1534 float phi = p.phi();
1535 return cdc.wire(r, phi);
1536}
1537
1538AList<TSegment>
1539TConformalFinder::pickUpSegmentsInConformal(float phi[12],
1540 float loadWidth,
1541 unsigned axialStereoSwitch) const {
1542 AList<TSegment> outList;
1543 HepPoint3D center(1., 1., 0.);
1544 center.setPhi(phi[0]);
1545
1546 //...Search for segments...
1547 for (unsigned sl = 0; sl < 2; sl++) {
1548 if (sl == 0) {
1549 if (! (axialStereoSwitch & 1)) continue;
1550 }
1551 else {
1552 if (! (axialStereoSwitch & 2)) continue;
1553 }
1554 for (unsigned i = 0; i < 6; i++) {
1555 const AList<TSegment> & list = _allUnused[sl][i];
1556 unsigned n = list.length();
1557 for (unsigned j = 0; j < n; j++) {
1558 TSegment & s = * list[j];
1559
1560 const HepPoint3D & p = s.position();
1561 if (center.dot(p) < 0.) continue;
1562
1563 bool found = false;
1564 const AList<TMLink> & inners = s.inners();
1565 unsigned m = inners.length();
1566 for (unsigned k = 0; k < m; k++) {
1567 float dPhi = phi[i * 2 + sl] -
1568 inners[k]->position().phi();
1569 dPhi = fabs(fmod(dPhi, PI2));
1570 if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
1571 if (dPhi < WIDTH[i * 2 + sl] * loadWidth) {
1572 outList.append(s);
1573 found = true;
1574 break;
1575 }
1576 }
1577
1578 if (found) continue;
1579 const AList<TMLink> & outers = s.outers();
1580 m = outers.length();
1581 for (unsigned k = 0; k < m; k++) {
1582 float dPhi = phi[i * 2 + sl + 1] -
1583 outers[k]->position().phi();
1584 dPhi = fabs(fmod(dPhi, PI2));
1585 if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
1586 if (dPhi < WIDTH[i * 2 + sl] * loadWidth) {
1587 outList.append(s);
1588 break;
1589 }
1590 }
1591 }
1592 }
1593 }
1594
1595 return outList;
1596}
1597
1598AList<TMLink>
1599TConformalFinder::pickUpLinksInConformal(float phi[12],
1600 float loadWidth,
1601 unsigned axialStereoSwitch) const {
1602 HepPoint3D center(1., 1., 0.);
1603 center.setPhi(phi[0]);
1604
1605 AList<TMLink> outList;
1606 unsigned nBad = _unused[2].length();
1607 for (unsigned i = 0; i < nBad; i++) {
1608 unsigned sl = _unused[2][i]->wire()->superLayerId();
1609 unsigned as = sl % 2;
1610 if (as == 0) {
1611 if (! (axialStereoSwitch & 1)) continue;
1612 }
1613 else {
1614 if (! (axialStereoSwitch & 2)) continue;
1615 }
1616
1617 const HepPoint3D & p = _unused[2][i]->position();
1618 if (center.dot(p) < 0.) continue;
1619
1620 float dPhi = phi[sl] - _unused[2][i]->position().phi();
1621 dPhi = fabs(fmod(dPhi, PI2));
1622 if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
1623 if (dPhi < WIDTH[sl] * loadWidth) {
1624 outList.append(_unused[2][i]);
1625 continue;
1626 }
1627 dPhi = phi[sl + 1] - _unused[2][i]->position().phi();
1628 dPhi = fabs(fmod(dPhi, PI2));
1629 if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
1630 if (dPhi < WIDTH[sl] * loadWidth)
1631 outList.append(_unused[2][i]);
1632 }
1633 return outList;
1634}
1635
1636int
1637TConformalFinder::crossPointsInConformal(const AList<TSegment> & inList,
1638 HepPoint3D points[12]) const {
1639 //...Parameters...
1640 static const float confRadius2[] = {4. / ( 8.3 * 8.3),
1641 4. / (16.9 * 16.9),
1642 4. / (21.7 * 21.7),
1643 4. / (31.3 * 31.3),
1644 4. / (36.1 * 36.1),
1645 4. / (44.1 * 44.1),
1646 4. / (50.5 * 50.5),
1647 4. / (58.5 * 58.5),
1648 4. / (64.9 * 64.9),
1649 4. / (72.9 * 72.9),
1650 4. / (79.3 * 79.3),
1651 4. / (87.4 * 87.4)};
1652
1653 //...Get conformal points from segments as seeds...
1654 AList<TMLink> forLine;
1655 HepPoint3D center;
1656 unsigned n = inList.length();
1657 for (unsigned i = 0; i < n; i++) {
1658 const HepPoint3D & p = inList[i]->position();
1659 forLine.append(new TMLink(NULL, NULL, p));
1660 center += p;
1661 }
1662 center *= 1. / float(n);
1663
1664 //...Make a line in the conformal plane...
1665 TMLine line(forLine);
1666 int err = line.fit();
1667 if (err) {
1668#ifdef TRKRECO_DEBUG
1669 std::cout << "crossPoints ... failed due to line fit failure" << std::endl;
1670#endif
1671 HepAListDeleteAll(forLine);
1672 return -1;
1673 }
1674
1675 //...Calculate points...
1676 for (unsigned i = 0; i < 12; i++) {
1677 HepPoint3D p[2];
1678 float c = - line.a() * line.b();
1679 float d = 1. + line.a() * line.a();
1680 float e = sqrt(confRadius2[i] * d - line.b() * line.b());
1681 p[0].setX((c + e) / d);
1682 p[0].setY(line.a() * p[0].x() + line.b());
1683 p[1].setX((c - e) / d);
1684 p[1].setY(line.a() * p[1].x() + line.b());
1685
1686 float mag0 = (center - p[0]).mag();
1687 float mag1 = (center - p[1]).mag();
1688
1689 if (mag0 < mag1) points[i] = p[0];
1690 else points[i] = p[1];
1691#ifdef TRKRECO_DEBUG_DETAIL
1692 std::cout << "0,1 = " << p[0] << ", " << p[1] << std::endl;
1693#endif
1694 }
1695
1696#ifdef TRKRECO_DEBUG_DETAIL
1697 std::cout << "... center is : " << center << std::endl;
1698 std::cout << "... cross points are : " << std::endl;
1699 for (unsigned i = 0; i < 12; i++) {
1700 std::cout << " " << i << " : " << points[i] << std::endl;
1701 }
1702#endif
1703#ifdef TRKRECO_WINDOW
1704 displayStatus("Conf::crossPoints ... line to salvage for conf plane");
1705 _rphiWindow.append(inList, leda_green);
1706 _rphiWindow.oneShot(line, leda_blue);
1707#endif
1708
1709 HepAListDeleteAll(forLine);
1710 return 0;
1711}
1712
1713AList<TSegment>
1714TConformalFinder::stereoSegments(const TTrack & t) const {
1715#ifdef TRKRECO_DEBUG
1716 std::cout << "... finding stereo segments" << std::endl;
1717#endif
1718
1719 const AList<TSegment> & bases = (const AList<TSegment> &) t.segments();
1720 AList<TSegment> seeds;
1721 unsigned n = bases.length();
1722 if (n == 0) return seeds;
1723
1724 //...Calculate cross points in the conformal plane...
1725 TPoint2D points[12];
1726 int err = crossPoints(t, points);
1727 //yuany
1728// for(unsigned i=0;i<12;i++)
1729// cout<<points[i]<<endl;
1730 if (err) {
1731#ifdef TRKRECO_DEBUG
1732 std::cout << "... stereo segment finding failed" << std::endl;
1733#endif
1734 return seeds;
1735 }
1736
1737 //...Pick up segments...
1738 // return pickUpSegments(points, float(_stereoLoadWidth), 2);
1739 AList<TSegment> list = pickUpSegments(points, float(_stereoLoadWidth), 2);
1740 //yuany
1741 //cout<<"list length "<<list.length()<<endl;
1742
1743 //...Save direction of a segment of axial layers...
1744 TPoint2D dir[5];
1745 for (unsigned i = 0; i < n; i++) {
1746 const TSegment & s = * bases[i];
1747 AList<TMLink> ll = s.links();
1748 unsigned sl = ll[0]->wire()->axialStereoLayerId()/4;
1749 dir[sl] = TPoint2D(s.direction());
1750// cout<<"in sl = "<<sl<<" dir = "<<dir[sl]<<endl;
1751 }
1752
1753 //...Cal. direction if empty...
1754 for (unsigned i = 0; i < 5; i++) {
1755 if (dir[i].mag() < .5) {
1756 unsigned j = i;
1757 while ((j < 5) && (dir[j].mag() < .5))
1758 ++j;
1759 if (j > 4) j = 4;
1760 if (dir[j].mag() < .5) {
1761 j = i;
1762 while ((j > 0) && (dir[j].mag() < .5))
1763 --j;
1764 }
1765 dir[i] = dir[j];
1766 }
1767 }
1768 //...Remove bad segments...
1769#ifdef TRKRECO_DEBUG_DETAIL
1770 std::cout << " ... direction :" << std::endl;
1771 for (unsigned i = 0; i < 6; i++)
1772 std::cout << " " << i << " : " << dir[i] << std::endl;
1773 std::cout << " ... direction cos :" << std::endl;
1774#endif
1775 AList<TSegment> badList;
1776 AList<TSegment> toBeChecked[6];
1777 unsigned nL = list.length();
1778 for (unsigned i = 0; i < 6; ++i)
1779 toBeChecked[i].removeAll();
1780
1781 for (unsigned j = 0; j < nL; ++j) {
1782 TSegment & s = * list[j];
1783 AList<TMLink> ll = s.links();
1784 unsigned sl = ll[0]->wire()->axialStereoLayerId()/4;
1785 toBeChecked[sl].append(s);
1786 }
1787
1788 for (unsigned i = 0; i < 6; ++i) {
1789 unsigned nToBeChecked = toBeChecked[i].length();
1790 if (nToBeChecked < 4) continue;
1791 unsigned nBadList = badList.length();
1792 for (unsigned j = 0; j < nToBeChecked; ++j) {
1793 if((badList.length() - nBadList) >= nToBeChecked - 4) break;
1794 TSegment & s = * toBeChecked[i][j];
1795 TPoint2D sDir = s.direction();
1796 unsigned sl = i;
1797 unsigned aSl;
1798
1799 if (sl < 2) aSl = 0;
1800 else if (sl < 4) aSl = 2;
1801 else aSl = 3;
1802 if (dir[aSl].dot(sDir) < 0.7) badList.append(s);
1803 if (sl == 2 || sl == 3 || sl == 4) {
1804 if (dir[aSl - 1].dot(sDir) < 0.6) badList.append(s);
1805 }
1806 else {
1807 if (dir[aSl + 1].dot(sDir) < 0.6) badList.append(s);
1808 }
1809 }
1810 toBeChecked[i].remove(badList);
1811 if (toBeChecked[i].length() > 4) {
1812 AList<TSegment> tmp; //with >= 4 links
1813 unsigned nSegs = toBeChecked[i].length();
1814 for (unsigned j = 0; j < nSegs; ++j) {
1815 TSegment & s = * toBeChecked[i][j];
1816 if (s.links().length() > 3) tmp.append(s);
1817 }
1818 toBeChecked[i].remove(tmp);
1819 for (unsigned j = 0; j < tmp.length(); ++j) {
1820 AList<TMLink> tmpLinks = tmp[j]->links();
1821 for (unsigned k = 0; k < toBeChecked[i].length(); ++k) {
1822 TSegment & ss = * toBeChecked[i][k];
1823 AList<TMLink> threeLinks = ss.links();
1824 for (unsigned kk = 0; kk < threeLinks.length(); ++kk) {
1825 TMLink & ll = * threeLinks[kk];
1826 if (tmpLinks.hasMember(ll)) {
1827 badList.append(ss);
1828 break;
1829 }
1830 }
1831 }
1832 }//j
1833 }
1834//jialk
1835 toBeChecked[i].remove(badList);
1836 if (toBeChecked[i].length() > 4) {//still with > 4 segs after above filter
1837 AList<TSegment> temp;
1838 unsigned nSegs = toBeChecked[i].length();
1839 for(unsigned k=0; k < nSegs; ++k){
1840 TSegment &s = *toBeChecked[i][k];
1841 temp.append(s);
1842 }
1843 //use bubble to get a sequence
1844 for(unsigned y=0; y<nSegs; y++){
1845 for(unsigned x=1; x<nSegs-y; x++){
1846 TSegment &s1 = *temp[x];
1847 TSegment &s2 = *temp[x-1];
1848 if(s1.links().length() < s2.links().length()){
1849 TSegment & s0 = *temp[x-1];
1850 *temp[x-1] = *temp[x];
1851 *temp[x] = s0;
1852 }
1853 }
1854 }
1855 for(unsigned j=4; j<nSegs; ++j){
1856 TSegment &s4 = *temp[j];
1857 badList.append(s4);
1858 }
1859 }
1860
1861
1862
1863 toBeChecked[i].remove(badList);
1864
1865// cout<<"**********the num of stereo segs after filter*********** "<<toBeChecked[i].length()<<endl;
1866 for(unsigned j = 0; j < toBeChecked[i].length(); ++j){
1867 TSegment &s5 = *toBeChecked[i][j];
1868// cout<<"num of hits"<<s5.links().length()<<endl;
1869 }
1870
1871 }
1872
1873/* for (unsigned i = 0; i < nL; i++) {
1874 TSegment & s = * list[i];
1875 AList<TMLink> ll = s.links();
1876 unsigned sl = ll[0]->wire()->axialStereoLayerId()/4;
1877 TPoint2D sDir = s.direction();
1878 unsigned aSl;
1879
1880 if (sl < 2) aSl = 0;
1881 else if (sl < 4) aSl = 2;
1882 else aSl = 3;
1883 if (dir[aSl].dot(sDir) < 0.7) badList.append(s);
1884 if (sl == 2 || sl == 3 || sl == 4) {
1885 if (dir[aSl - 1].dot(sDir) < 0.6) badList.append(s);
1886 }
1887 else {
1888 if (dir[aSl + 1].dot(sDir) < 0.6) badList.append(s);
1889 }
1890#ifdef TRKRECO_DEBUG_DETAIL
1891 std::cout << " " << dir[sl].dot(sDir) << ", ";
1892 std::cout << dir[sl + 1].dot(sDir) << " : ";
1893 s.dump("hits sort");
1894#endif
1895 }
1896*/
1897 //...Width cut...
1898// for (unsigned i = 0; i < nL; i++) {
1899// TSegment & s = * list[i];
1900// unsigned width = s.width();
1901// if (width > 2) badList.append(s);
1902// }
1903
1904// cout<<"stereoSegment: "<<list.length()<<" badList: "<<badList.length()<<endl;
1905 list.remove(badList);
1906// cout<<"after remove bads: "<<list.length()<<endl;
1907
1908#ifdef TRKRECO_DEBUG
1909 std::cout << " ... bad segments :" << std::endl;
1910 for (unsigned i = 0; i < badList.length(); i++)
1911 badList[i]->dump("hits sort", " ");
1912#endif
1913
1914 return list;
1915}
1916
1917HepVector3D
1918TConformalFinder::direction(const TSegment & segment) const {
1919 AList<TMLink> forLine;
1920 const TSegment * s = & segment;
1921 while (1) {
1922 if (s->outerLinks().length() != 1) break;
1923 const TSegment * o = s->outerLinks()[0];
1924 const HepPoint3D & p = o->position();
1925 forLine.append(new TMLink(NULL, NULL, p));
1926 s = o;
1927 }
1928
1929 if (forLine.length() == 0)
1930 return segment.direction();
1931 else if (forLine.length() < 2) {
1932 HepAListDeleteAll(forLine);
1933 return segment.direction();
1934 }
1935
1936 //...Make a line in the conformal plane...
1937 TMLine line(forLine);
1938 int err = line.fit();
1939 if (err) {
1940#ifdef TRKRECO_DEBUG
1941 std::cout << "direction ... failed due to line fit failure" << std::endl;
1942#endif
1943 HepAListDeleteAll(forLine);
1944 return segment.direction();
1945 }
1946
1947 HepAListDeleteAll(forLine);
1948 Vector3 tmp(-1., -(line.a() + line.b()));
1949 tmp.unit();
1950 return HepVector3D(tmp);
1951}
1952
1953
1954void
1955TConformalFinder::salvage(TTrack & track,
1956 unsigned asSwitch,
1957 const AList<TSegment> & bads) const {
1958#ifdef TRKRECO_DEBUG_DETAIL
1959 std::cout << "... salvaging in real plane" << std::endl;
1960#endif
1961
1962 //...Calculate cross points...
1963 TPoint2D points[12];
1964 int err = crossPoints(track, points);
1965 if (err) {
1966#ifdef TRKRECO_DEBUG_DETAIL
1967 std::cout << "... salvaging failed in calculation of intersection"
1968 << std::endl;
1969#endif
1970 return;
1971 }
1972
1973 //...Pick up segments...
1974 AList<TSegment> toBeChecked = pickUpSegments(points,
1975 float(_salvageLoadWidth),
1976 asSwitch);
1977 toBeChecked.remove(bads);
1978 toBeChecked.remove(track.segments());
1979 toBeChecked = trackSide(track, toBeChecked);
1980
1981 //...Pick up links...
1982 AList<TMLink> links;
1983 unsigned nB = toBeChecked.length();
1984 for (unsigned i = 0; i < nB; i++) {
1985 const AList<TMLink> & tmp = toBeChecked[i]->links();
1986 unsigned nL = tmp.length();
1987 for (unsigned j = 0; j < nL; j++) {
1988 if (tmp[j]->hit()->track()) continue;
1989 links.append(tmp[j]);
1990 }
1991 }
1992
1993 //...Search in bad hits...
1994 AList<TMLink> all = pickUpLinks(points, float(_salvageLoadWidth), asSwitch);
1995#ifdef TRKRECO_DEBUG
1996 std::cout<<" pickUpLinks length = "<<all.length()<<std::endl;
1997#endif
1998
1999 all = trackSide(track, all);
2000 all.remove(track.links());
2001#ifdef TRKRECO_DEBUG
2002 std::cout<<" after reomove pickUpLinks length = "<<all.length()<<std::endl;
2003#endif
2004 links.append(all);
2005
2006#ifdef TRKRECO_WINDOW
2007 AList<TMLink> tmp = links;
2008#endif
2009
2010 //...Ask builder...
2011 _builder.salvage(track, links);
2012
2013 //...Check used segments...
2014 const AList<TMLink> & usedLinks = track.links();
2015 for (unsigned i = 0; i < nB; i++) {
2016 TSegment & segment = * toBeChecked[i];
2017 AList<TMLink> used = Links(segment, track);
2018 if (used.length() == 0) continue;
2019
2020 track.segments().append(segment);
2021 segment.tracks().append(track);
2022 }
2023
2024#ifdef TRKRECO_WINDOW
2025 displayStatus("Conf::salvage ... results");
2026 TTrackBase y(tmp);
2027 _rphiWindow.append(y, leda_red);
2028 _rphiWindow.append(toBeChecked, leda_red);
2029 _rphiWindow.oneShot(track, leda_green);
2030#endif
2031
2032 //...Termination...
2033 return;
2034}
2035
2036int
2037TConformalFinder::crossPoints(const TTrack & track,
2038 TPoint2D points[12]) const {
2039
2040#ifdef TRKRECO_DEBUG_DETAIL
2041 std::cout << " cross points are :" << std::endl;
2042#endif
2043
2044 /*
2045 //...Parameters...
2046 static const float R[] = { 8.3,
2047 16.9,
2048 21.7,
2049 31.3,
2050 36.1,
2051 44.1,
2052 50.5,
2053 58.5,
2054 64.9,
2055 72.9,
2056 79.3,
2057 87.4};
2058 static const float R2[] = { 8.3 * 8.3,
2059 16.9 * 16.9,
2060 21.7 * 21.7,
2061 31.3 * 31.3,
2062 36.1 * 36.1,
2063 44.1 * 44.1,
2064 50.5 * 50.5,
2065 58.5 * 58.5,
2066 64.9 * 64.9,
2067 72.9 * 72.9,
2068 79.3 * 79.3,
2069 87.4 * 87.4};
2070 */
2071 static const float R[] = { 7.3,
2072 12.21,
2073 18.9,
2074 25.375,
2075 31.86,
2076 39.16,
2077 45.685,
2078 52.215,
2079 58.665,
2080 65.91,
2081 72.37,
2082 80.};
2083 static const float R2[] = { R[0]*R[0],
2084 R[1]*R[1],
2085 R[2]*R[2],
2086 R[3]*R[3],
2087 R[4]*R[4],
2088 R[5]*R[5],
2089 R[6]*R[6],
2090 R[7]*R[7],
2091 R[8]*R[8],
2092 R[9]*R[9],
2093 R[10]*R[10],
2094 R[11]*R[11]};
2095
2096
2097 static const TPoint2D o(0., 0.);
2098 TPoint2D c = track.center();
2099 TPoint2D co = - c;
2100 double r = fabs(track.radius());
2101 double r2 = r * r;
2102 double d2 = c.mag2();
2103 double d = sqrt(d2);
2104 double sl = - c.x() / c.y();
2105 //...Calculate points...
2106 unsigned nOk = 0;
2107 for (unsigned i = 0; i < 12; i++) {
2108 double minR = r < R[i] ? r : R[i];
2109 double maxR = r < R[i] ? R[i] : r;
2110
2111#ifdef TRKRECO_DEBUG_DETAIL
2112 std::cout << " r,R ... " << minR << ", " << maxR << ", " << d << std::endl;
2113#endif
2114
2115 if ((r + R[i] < d) || (minR + d < maxR)) { //no intersection of two circles
2116 points[i] = o;
2117 continue;
2118 }
2119 ++nOk;
2120 double a = R2[i] + d2 - r2;
2121 double s = sqrt(4. * R2[i] * d2 - a * a);
2122 double q = 0.5 * a / c.y();
2123 points[i].x(0.5 * (c.x() * a + c.y() * s) / d2);
2124 points[i].y(q + sl * points[i].x());
2125 const double Bz = -1000*(m_pmgnIMF->getReferField());
2126 if (co.cross(points[i] - c) * track.charge() * Bz > 0.) {
2127 points[i].x(0.5 * (c.x() * a - c.y() * s) / d2);
2128 points[i].y(q + sl * points[i].x());
2129 }
2130#ifdef TRKRECO_DEBUG_DETAIL
2131 std::cout << " " << i << " : " << points[i] << std::endl;
2132 std::cout << " chg=" << track.charge()<<std::endl;
2133 std::cout << ", c=" << c << ", co=" <<std::endl;
2134 std::cout << ", " << co.cross(points[i] - c) * track.charge() << std::endl;
2135 std::cout << " " << 0.5 * (c.x() * a + c.y() * s) / d2;
2136 std::cout << ", " << q + sl * (0.5 * (c.x() * a + c.y() * s) / d2);
2137 std::cout << " " << 0.5 * (c.x() * a - c.y() * s) / d2;
2138 std::cout << ", " << q + sl * (0.5 * (c.x() * a - c.y() * s) / d2);
2139 std::cout << std::endl;
2140#endif
2141
2142 }
2143 if (nOk) return 0;
2144 else return -1;
2145}
2146
2147AList<TSegment>
2148TConformalFinder::pickUpSegments(const TPoint2D x[12],
2149 float loadWidth,
2150 unsigned axialStereoSwitch) const {
2151 static const TPoint2D O(0., 0.);
2152 AList<TSegment> outList;
2153
2154 //yuany
2155 //...Search for segments...
2156 for (unsigned sl = 0; sl < 2; sl++) {
2157 unsigned nMax = 5;
2158 if (sl == 0) {
2159 if (! (axialStereoSwitch & 1)) continue;
2160 }
2161 else {
2162 nMax = 6;
2163 if (! (axialStereoSwitch & 2)) continue;
2164 }
2165 for (unsigned i = 0; i < nMax; i++) {
2166
2167 unsigned layerId;
2168 switch(sl){
2169 case 0:
2170 switch(i){
2171 case 0:layerId = 2;break;
2172 case 1:layerId = 3;break;
2173 case 2:layerId = 4;break;
2174 case 3:layerId = 9;break;
2175 case 4:layerId = 10;break;
2176 }
2177 case 1:
2178 switch(i){
2179 case 0:layerId = 0;break;
2180 case 1:layerId = 1;break;
2181 case 2:layerId = 5;break;
2182 case 3:layerId = 6;break;
2183 case 4:layerId = 7;break;
2184 case 5:layerId = 8;break;
2185 }
2186 }
2187 // cout<<"sl "<<sl<<" i "<<i<<" layerId "<<layerId<<endl;
2188
2189 if (x[layerId] == O) continue;
2190 float a = WIDTH[layerId] * loadWidth;
2191 float phi0 = x[layerId].phi();
2192 float phi1 = x[layerId+1].phi();//Liuqg
2193 float phiC = (phi0 + phi1) / 2.;
2194 float phiD = fabs(phi0 - phiC);
2195 if (phiD > M_PI_2) {
2196 phiC += M_PI;
2197 phiD = M_PI - phiD;
2198 }
2199 phiD += a;
2200
2201 const AList<TSegment> & list = _allUnused[sl][i];
2202 unsigned n = list.length();
2203#ifdef TRKRECO_DEBUG_DETAIL
2204 std::cout << " pickUpSegments ... phi0,1,C,D=" << phi0 << ",";
2205 std::cout << phi1 << "," << phiC << "," << phiD <<" nhits "<<n<< std::endl;
2206#endif
2207 for (unsigned j = 0; j < n; j++) {
2208 TSegment & s = * list[j];
2209
2210#ifdef TRKRECO_DEBUG_DETAIL
2211 s.dump("hits sort", " ");
2212#endif
2213
2214 bool found = false;
2215 const AList<TMLink> & inners = s.inners();
2216 unsigned m = inners.length();
2217 for (unsigned k = 0; k < m; k++) {
2218 float phi = inners[k]->position().phi();
2219 if (phi < 0.) phi += PI2;
2220 float dPhi = fabs(phi - phiC);
2221 if (dPhi > M_PI) dPhi = PI2 - dPhi;
2222#ifdef TRKRECO_DEBUG_DETAIL
2223 std::cout << " " << inners[k]->wire()->name();
2224 std::cout << ", phi,dPhi=" << phi << "," << dPhi << std::endl;
2225#endif
2226 if (dPhi < phiD) {
2227 outList.append(s);
2228 found = true;
2229 break;
2230 }
2231 }
2232
2233 if (found) continue;
2234 const AList<TMLink> & outers = s.outers();
2235 m = outers.length();
2236 for (unsigned k = 0; k < m; k++) {
2237 float phi = outers[k]->position().phi();
2238 if (phi < 0.) phi += PI2;
2239 float dPhi = fabs(phi - phiC);
2240 if (dPhi > M_PI) dPhi = PI2 - dPhi;
2241#ifdef TRKRECO_DEBUG_DETAIL
2242 std::cout << " " << outers[k]->wire()->name();
2243 std::cout << ", phi,dPhi=" << phi << "," << dPhi << std::endl;
2244#endif
2245 if (dPhi < phiD) {
2246 outList.append(s);
2247 found = true;
2248 break;
2249 }
2250 }
2251 }
2252 }
2253 }
2254 return outList;
2255}
2256
2257AList<TMLink>
2258TConformalFinder::pickUpLinks(const TPoint2D x[12],
2259 float loadWidth,
2260 unsigned axialStereoSwitch) const {
2261
2262 static const TPoint2D O(0., 0.);
2263 AList<TMLink> outList;
2264
2265 unsigned nBad = _unused[2].length();
2266 for (unsigned i = 0; i < nBad; i++) {
2267 unsigned sl = _unused[2][i]->wire()->superLayerId();
2268 unsigned as = 1;
2269 if (_unused[2][i]->wire()->axial()) as = 0;
2270 if (as == 0) {
2271 if (! (axialStereoSwitch & 1)) continue;
2272 }
2273 else {
2274 if (! (axialStereoSwitch & 2)) continue;
2275 }
2276
2277 float a = WIDTH[sl] * loadWidth;
2278 float phi0 = x[sl].phi();
2279 float phi1 = x[sl + 1].phi();
2280 float phi = _unused[2][i]->position().phi();
2281
2282 if (phi < 0.) phi += PI2;
2283 if (phi1 < phi0) {
2284 phi1 = phi0;
2285 phi0 = x[sl + 1].phi();
2286 }
2287 float dPhi = phi1 - phi0;
2288 if (dPhi < M_PI) {
2289 phi0 -= a;
2290 phi1 += a;
2291 if (phi > phi0 && phi < phi1) outList.append(_unused[2][i]);
2292 }
2293 else {
2294 phi0 += a;
2295 phi1 -= a;
2296 if (phi < phi0 || phi > phi1) outList.append(_unused[2][i]);
2297 }
2298#ifdef TRKRECO_DEBUG_DETAIL
2299 std::cout << "TConformalFinder:: pickUpLinks " << std::endl;
2300 std::cout << "a "<<a <<" phi0 "<<phi0<<" phi1 "<<phi1<<" phi "<< phi<<" dPhi "<<dPhi << std::endl;
2301#endif
2302 }
2303 return outList;
2304}
2305
2306void
2307TConformalFinder::slowFinding2D(unsigned level) {
2308
2309#ifdef TRKRECO_DEBUG_DETAIL
2310 std::cout << "... slow finding (2D) : level = " << level << std::endl;
2311#endif
2312#ifdef TRKRECO_DEBUG
2313 _stage = ConformalSlow2D;
2314#endif
2315
2316 AList<TSegment> * segments = & _allUnused[0][0];
2317
2318 unsigned id = _2DTracks.length() + _3DTracks.length();
2319 unsigned seedLayer = 5;
2320 while (seedLayer--) {
2321
2322 //...Seed loop...
2323 AList<TSegment> tmp = segments[seedLayer];
2324 unsigned n = tmp.length();
2325 for (unsigned i = 0; i < n; i++) {
2326 AList<TTrack> trackCandidates;
2327 TSegment & current = * tmp[i];
2328 if (current.links().length() < 3) continue;
2329 if (current.innerLinks().length() != 1) continue;
2330 if (current.innerLinks()[0]->links().length() < 3) continue;
2331 AList<TSegment> seeds;
2332 seeds.append(current);
2333 seeds.append(current.innerLinks()[0]);
2334
2335 std::string name = "s" + itostring(id);
2336
2337#ifdef TRKRECO_DEBUG
2338 std::cout << "==> slow building : " << name << std::endl;
2339#endif
2340
2341 //...Try to build a track...
2342 for (int k = 0; k < seeds.length(); ++k) //tsf
2343 seeds[k]->solveLR();
2344
2345 TTrack * t = expand(seeds);
2346 if (t) {
2347 AList<TSegment> bads;
2348 t->name(name);
2349 salvage(* t, 1, bads);
2350 if (! trackQuality(* t)) {
2351 deleteTrack(* t);
2352 continue;
2353 }
2354 t->finder(TrackSlowFinder);
2355 t->quality(TrackQuality2D);
2356 t->fitting(TrackFitGlobal);
2357 trackCandidates.append(t);
2358
2359#ifdef TRKRECO_DEBUG
2360 std::cout << "... 2D finished : " << t->name() << std::endl;
2361 t->dump("hits sort flag pull", " ");
2362#endif
2363#ifdef TRKRECO_WINDOW
2364 displayStatus("Conf::expand ... finished");
2365 _rphiWindow.oneShot(* t, leda_green);
2366#endif
2367 removeUsedSegments(trackCandidates);
2368 _2DTracks.append(trackCandidates);
2369 id = _2DTracks.length() + _3DTracks.length();
2370 }
2371 }
2372
2373 //...Resolve multi-track candidates...
2374 // resolveSegments(trackCandidates);
2375
2376 //...Remove used segments...
2377 }
2378}
2379
2380TTrack *
2381TConformalFinder::expand(AList<TSegment> & seeds) const {
2382#ifdef TRKRECO_WINDOW
2383 displayStatus("Conf::expand ... seeds");
2384 _rphiWindow.oneShot(seeds, leda_green);
2385#endif
2386#ifdef TRKRECO_DEBUG
2387 std::cout << "... expand : # of seeds = " << seeds.length() << std::endl;
2388#endif
2389 TTrack * t = NULL;
2390 if (seeds.length() > 2) {
2391 TTrack * t = _builder.buildRphi(seeds);
2392 if (t) {
2393 if (! trackQuality(* t)) {
2394 deleteTrack(* t);
2395 t = NULL;
2396 }
2397 }
2398 }
2399 if (t == NULL) {
2400 AList<TMLink> links = Links(seeds);
2401 TCircle c(links);
2402 c.fit();
2403 t = new TTrack(c);
2404 t->fit();
2405 t->segments().append(seeds);
2406 }
2407 return expand(* t);
2408}
2409
2410TTrack *
2411TConformalFinder::expand(TTrack & ti) const {
2412#ifdef TRKRECO_WINDOW
2413 displayStatus("Conf::expand ... seed track");
2414 _rphiWindow.oneShot(ti, leda_green);
2415#endif
2416#ifdef TRKRECO_DEBUG
2417 std::cout << "... expand : by a track" << std::endl;
2418#endif
2419
2420 TTrack * t = & ti;
2421 AList<TSegment> seeds = t->segments();
2422 unsigned nSegments = seeds.length();
2423 unsigned nCores = t->nCores();
2424 unsigned nTrial = 0;
2425 while (1) {
2426
2427 //...Target superlayer...
2428 unsigned sl = SuperLayer(t->cores());
2429 unsigned inner;
2430 unsigned outer;
2431 targetSuperLayer(sl, inner, outer);
2432//cout<<"sl: "<<sl<<" inner: "<<inner<<" outer: "<<outer<<endl;
2433 unsigned target = inner;
2434 if (target == 11)
2435 target = outer;
2436 if (target == 11)
2437 break;
2438
2439 //...Get target segments...
2440 AList<TSegment> segments = targetSegments(* t, target);
2441 if (segments.length() == 0) {
2442 unsigned fl[5] = {2, 3, 4, 9, 10};
2443 unsigned inner2 = 0;
2444 unsigned outer2 = 5;
2445 for (int i = 0; i < 5; i++) {
2446 if (fl[i] == inner) inner2 = i;
2447 if (fl[i] == outer) outer2 = i;
2448 }
2449 if (inner > 2&& inner2>0) target = fl[inner2 - 1];
2450 if (target == 9) {
2451 target = fl[outer2 + 1];
2452 }
2453 if (target == 13) break;
2454 segments = targetSegments(* t, target);
2455 }
2456#ifdef TRKRECO_DEBUG
2457 std::cout << "... expand : segments to be checked = ";
2458 std::cout << segments.length() << std::endl;
2459#endif
2460#ifdef TRKRECO_WINDOW
2461 displayStatus("Conf::expand ... candidate segments");
2462 _rphiWindow.append(segments, leda_blue);
2463 _rphiWindow.oneShot(* t, leda_green);
2464#endif
2465
2466 //...Create candidate tracks...
2467 unsigned n = segments.length();
2468 AList<TTrack> tracks;
2469 for (unsigned i = 0; i < n; i++) {
2470 segments[i]->solveLR();
2471 AList<TSegment> seed0 = seeds;
2472 seed0.append(segments[i]);
2473 TTrack * t0 = _builder.buildRphi(seed0);
2474 if (! t0) continue;
2475 if ((t0->segments().length() <= nSegments) ||
2476 (t0->nCores() <= nCores) ||
2477 (SuperLayer(t0->cores()) <= sl)) {
2478 deleteTrack(* t0);
2479 continue;
2480#ifdef TRKRECO_DEBUG
2481 std::cout << "... expand : candidate deleted" << std::endl;
2482#endif
2483 }
2484 tracks.append(t0);
2485 }
2486
2487#ifdef TRKRECO_DEBUG
2488 std::cout << "... expand : target superlayer = " << target;
2489 std::cout << ", # of track candidates = " << tracks.length()
2490 << std::endl;
2491#endif
2492
2493 //...Expand by links...
2494 if (tracks.length() == 0) {
2495 AList<TMLink> links = targetLinks(* t, target);
2496 _builder.salvage(* t, links);
2497 unsigned newSl = SuperLayer(t->cores());
2498 if (newSl & (1 << (target))) {
2499#ifdef TRKRECO_DEBUG
2500 std::cout << "... expand : links to be appended = ";
2501 std::cout << links.length() << std::endl;
2502#endif
2503#ifdef TRKRECO_WINDOW
2504 displayStatus("Conf::expand ... candidate links");
2505 _rphiWindow.append(links, leda_blue);
2506 _rphiWindow.oneShot(* t, leda_green);
2507#endif
2508 TTrack * t0 = new TTrack(* t);
2509 if (! t0) break;
2510 if ((t0->nCores() <= nCores) ||
2511 (SuperLayer(t0->cores()) <= sl)) {
2512 deleteTrack(* t0);
2513 break;
2514#ifdef TRKRECO_DEBUG
2515 std::cout << "... expand : candidate deleted" << std::endl;
2516#endif
2517 }
2518 tracks.append(t0);
2519 }
2520 else
2521 break;
2522 }
2523
2524#ifdef TRKRECO_DEBUG
2525 std::cout << "... expand : # of track candidates = " << tracks.length()
2526 << std::endl;
2527#endif
2528
2529 //...Select best...
2530 unsigned nTracks = tracks.length();
2531 if (nTracks == 0) break;
2532 TTrack * tNew = tracks[0];
2533 unsigned nBestCores = tNew->nCores();
2534 for (unsigned i = 1; i < nTracks; i++) {
2535 if (tracks[i]->nCores() > nBestCores) {
2536 deleteTrack(* tNew);
2537 tNew = tracks[i];
2538 nBestCores = tNew->nCores();
2539 }
2540 else {
2541 deleteTrack(* tracks[i]);
2542 }
2543 }
2544 nSegments = tNew->segments().length();
2545 nCores = nBestCores;
2546 seeds = tNew->segments();
2547
2548 //...Quality check...
2549 if (! trackQuality(* tNew)) {
2550 deleteTrack(* tNew);
2551 break;
2552 }
2553 deleteTrack(* t);
2554 t = tNew;
2555 }
2556
2557#ifdef TRKRECO_DEBUG
2558 if (t == & ti) std::cout << "... expand : failed" << std::endl;
2559 else std::cout << "... expand : track expanded" << std::endl;
2560#endif
2561#ifdef TRKRECO_WINDOW
2562 displayStatus("Conf::expand ... finished");
2563 _rphiWindow.oneShot(* t, leda_green);
2564#endif
2565
2566 return t;
2567}
2568
2569void
2570TConformalFinder::targetSuperLayer(unsigned sl,
2571 unsigned & inner,
2572 unsigned & outer) const {
2573 inner = 11;
2574 outer = 11;
2575 bool innerFound = false;
2576 bool outerFound = false;
2577 for (unsigned i = 0; i < 5; i++) { //changed to BESIII,Liuqg
2578 unsigned fl[5] ={2, 3, 4, 9, 10};
2579 if (! innerFound) {
2580 if (sl & (1 << (fl[i]))) innerFound = true;
2581 else inner = fl[i];
2582 }
2583 if (! outerFound) {
2584 if (sl & (1 << (fl[4 - i]))) outerFound = true;
2585 else outer = fl[4 - i];
2586 }
2587 }
2588}
2589
2590AList<TSegment>
2591TConformalFinder::targetSegments(const TTrack & t, unsigned sl) const {
2592 AList<TSegment> targets;
2593
2594 //...Calculate cross points...
2595 TPoint2D points[12];
2596 int err = crossPoints(t, points);
2597 if (err) {
2598#ifdef TRKRECO_DEBUG
2599 std::cout << " ... no target found" << std::endl;
2600#endif
2601 return targets;
2602 }
2603
2604 //...Pick up segments...
2605 AList<TSegment> toBeChecked = pickUpSegments(points, 3, 1);
2606 unsigned n = toBeChecked.length();
2607 for (unsigned i = 0; i < n; i++) {
2608 if (toBeChecked[i]->superLayerId() == sl)
2609 targets.append(toBeChecked[i]);
2610 }
2611
2612 return targets;
2613}
2614
2615AList<TMLink>
2616TConformalFinder::targetLinks(const TTrack & t, unsigned sl) const {
2617 AList<TMLink> targets;
2618
2619 //...Calculate cross points...
2620 TPoint2D points[12];
2621 int err = crossPoints(t, points);
2622 if (err) {
2623#ifdef TRKRECO_DEBUG
2624 std::cout << " ... no target found" << std::endl;
2625#endif
2626 return targets;
2627 }
2628
2629 //...Pick up segments...
2630 AList<TMLink> toBeChecked = pickUpLinks(points, 3, 1);
2631 unsigned n = toBeChecked.length();
2632 for (unsigned i = 0; i < n; i++) {
2633 if (toBeChecked[i]->wire()->superLayerId() == sl)
2634 targets.append(toBeChecked[i]);
2635 }
2636
2637 return targets;
2638}
2639
2640AList<TSegment>
2641TConformalFinder::refineSegments(const TTrack & t) const {
2642 const AList<TSegment> & original = t.segments();
2643 AList<TSegment> outList;
2644 unsigned n = original.length();
2645 for (unsigned i = 0; i < n; i++) {
2646 AList<TMLink> tmp = Links(* original[i], t);
2647 if (tmp.length() >= _minNLinksForSegmentInRefine)
2648 outList.append(original[i]);
2649 }
2650 unsigned sl = SuperLayer(outList);
2651 unsigned nCMax = 0;
2652 unsigned nStart = 0;
2653 unsigned nC = 0;
2654 unsigned nS = 0;
2655 unsigned fl[5] = {2, 3, 4, 9, 10};
2656 for (unsigned i = 0; i < 5; i++) {
2657 if (sl & (1 << fl[i])) {
2658 ++nC;
2659 if (nC == 1) nS = i;
2660 if (nC > nCMax) {
2661 nCMax = nC;
2662 nStart = nS;
2663 }
2664 }
2665 else {
2666 nC = 0;
2667 }
2668 }
2669// nStart *= 2;
2670// nCMax *= 2;
2671
2672 outList.removeAll();
2673 if (nCMax >= _minNSegments) {
2674 for (unsigned i = 0; i < n; i++) {
2675 unsigned id = original[i]->superLayerId();
2676 if ((id >= fl[nStart]) && (id < fl[nStart + nCMax]))
2677 outList.append(original[i]);
2678 }
2679 }
2680
2681#ifdef TRKRECO_DEBUG
2682 std::cout << "... refine segments : orignal sl = ";
2683 bitDisplay(sl, 11, 0);
2684 std::cout << ", output sl = ";
2685 bitDisplay(SuperLayer(outList), 11, 0);
2686 std::cout << std::endl;
2687#endif
2688
2689 return outList;
2690}
2691
2692bool
2693TConformalFinder::trackQuality(const TTrack & t) const {
2694 unsigned n1 = NSuperLayers(t.links());
2695 unsigned n3 = NSuperLayers(t.links(), 2);
2696
2697 //...check # radius...Liuqg
2698 if(fabs(t.radius())<15.) return false;
2699
2700#ifdef TRKRECO_WINDOW
2701 std::cout << "... trackQuality : n1,n3,nMissing=" << n1 << "," << n3;
2702 std::cout << "," << NMissingAxialSuperLayers(t.links()) << std::endl;
2703#endif
2704#ifdef TRKRECO_WINDOW
2705 displayStatus("trackQuality");
2706 if ((n3 < 2) || (NMissingAxialSuperLayers(t.links()) > 1))
2707 _rphiWindow.oneShot(t, leda_red);
2708 else
2709 _rphiWindow.oneShot(t, leda_green);
2710#endif
2711 if (n3 < 2) return false;
2712// if (NMissingAxialSuperLayers(t.links()) > 1) return false;
2713
2714 return true;
2715}
2716
2717void
2718TConformalFinder::refineLinks(TTrack & t, unsigned minN) const {
2719 const AList<TMLink> & links = t.links();
2720 AList<TMLink> sl[11];
2721 unsigned n = links.length();
2722 for (unsigned i = 0; i < n; i++)
2723 sl[links[i]->wire()->superLayerId()].append(links[i]);
2724 AList<TMLink> toBeRemoved;
2725 for (unsigned i = 0; i < 11; i++)
2726 if (sl[i].length() < minN)
2727 toBeRemoved.append(sl[i]);
2728 t.remove(toBeRemoved);
2729#ifdef TRKRECO_DEBUG
2730 std::cout << "... refining links : removed links = " << toBeRemoved.length();
2731 std::cout << std::endl;
2732#endif
2733}
2734
2735AList<TMLink>
2736TConformalFinder::trackSide(const TTrack & t, const AList<TMLink> & a) const {
2737 static const TPoint2D o(0., 0.);
2738 TPoint2D c = t.center();
2739 TPoint2D co = - c;
2740
2741 AList<TMLink> tmp;
2742 unsigned n = a.length();
2743 for (unsigned i = 0; i < n; i++) {
2744 TPoint2D x = a[i]->wire()->xyPosition();
2745 const double Bz = -1000*(m_pmgnIMF->getReferField());
2746 if (co.cross(x - c) * t.charge() * Bz < 0.)
2747 tmp.append(a[i]);
2748 }
2749 return tmp;
2750}
2751
2752AList<TSegment>
2753TConformalFinder::trackSide(const TTrack & t,
2754 const AList<TSegment> & a) const {
2755 static const TPoint2D o(0., 0.);
2756 TPoint2D c = t.center();
2757 TPoint2D co = - c;
2758
2759 AList<TSegment> tmp;
2760 unsigned n = a.length();
2761 for (unsigned i = 0; i < n; i++) {
2762 const AList<TMLink> & b = a[i]->links();
2763 bool bad = false;
2764 unsigned nB = b.length();
2765 for (unsigned j = 0; j < nB; j++) {
2766 TPoint2D x = b[j]->wire()->xyPosition();
2767 //yzhang 2012-05-03
2768 const double Bz = -1000*(m_pmgnIMF->getReferField());
2769 if (co.cross(x - c) * t.charge()* Bz > 0.) {
2770 bad = true;
2771 break;
2772 }
2773 }
2774 if (bad) continue;
2775 tmp.append(a[i]);
2776 }
2777 return tmp;
2778}
2779
2780void
2781TConformalFinder::resolveHits(AList<TTrack> & tracks) const {
2782 unsigned nTracks = tracks.length();
2783 if (nTracks < 2) return;
2784
2785 for (unsigned i = 0; i < nTracks - 1; i++) {
2786 TTrack & t0 = * tracks[i];
2787 const AList<TMLink> & links0 = t0.links();
2788 unsigned n0 = links0.length();
2789 AList<TMLink> remove0;
2790 for (unsigned j = i + 1; j < nTracks; j++) {
2791 TTrack & t1 = * tracks[j];
2792 const AList<TMLink> & links1 = t1.links();
2793 unsigned n1 = links1.length();
2794 AList<TMLink> remove1;
2795
2796 //...Check overlap...
2797 for (unsigned k = 0; k < n0; k++) {
2798 TMLink & l = * links0[k];
2799
2800 //...Decide which is better...
2801 if (links1.hasMember(l)) {
2802 TMLink l0(NULL, l.hit());
2803 TMLink l1(NULL, l.hit());
2804 t0.approach(l0);
2805 t1.approach(l1);
2806 if (l0.pull() > l1.pull())
2807 remove0.append(l);
2808 else
2809 remove1.append(l);
2810 }
2811 }
2812
2813 if (remove1.length()) {
2814 t1.remove(remove1);
2815 t1.fit();
2816 }
2817 }
2818
2819 if (remove0.length()) {
2820 t0.remove(remove0);
2821 t0.fit();
2822 }
2823 }
2824}
2825
2826AList<TSegment>
2827TConformalFinder::stereoSegmentsFromBadHits(const TTrack & t) const {
2828
2829 AList<TSegment> output;
2830 for (unsigned i = 0; i < 6; i++)
2831 output.append(new TSegment());
2832
2833 //...Calculate cross points...
2834 TPoint2D points[12];
2835 int err = crossPoints(t, points);
2836 if (err) {
2837#ifdef TRKRECO_DEBUG_DETAIL
2838 std::cout << "... conf::stereoSegBadHits : "
2839 << "error in calculation of intersection" << std::endl;
2840#endif
2841 return output;
2842 }
2843
2844 static const TPoint2D O(0., 0.);
2845 unsigned nBad = _unused[2].length();
2846 for (unsigned i = 0; i < nBad; i++) {
2847 unsigned sl = _unused[2][i]->wire()->superLayerId();
2848 unsigned asl = _unused[2][i]->wire()->axialStereoLayerId()/4;
2849 unsigned as;
2850 if (_unused[2][i]->wire()->axial())
2851 as = 0;
2852 else as = 1;
2853// unsigned as = sl % 2;
2854 if (as == 0) continue;
2855
2856 float a = WIDTH[sl] * _salvageLoadWidth;
2857 float phi0 = points[sl].phi();
2858 float phi1 = points[sl + 1].phi();
2859 float phi = _unused[2][i]->position().phi();
2860 if (phi < 0.) phi += PI2;
2861 if (phi1 < phi0) {
2862 phi1 = phi0;
2863 phi0 = points[sl + 1].phi();
2864 }
2865 float dPhi = phi1 - phi0;
2866 if (dPhi < M_PI) {
2867 phi0 -= a;
2868 phi1 += a;
2869 if (phi > phi0 && phi < phi1)
2870// output[sl / 2]->append(* _unused[2][i]);
2871 output[asl / 4]->append(* _unused[2][i]);
2872 }
2873 else {
2874 phi0 += a;
2875 phi1 -= a;
2876 if (phi < phi0 || phi > phi1)
2877 output[asl / 4]->append(* _unused[2][i]);
2878 }
2879 }
2880
2881 return output;
2882}
2883
2884void
2885TConformalFinder::findSegmentsPerfect(void) {
2886
2887 //...Create lists of links for each super layer...
2888 AList<TMLink> links[11];
2889 unsigned nHits = _hits[2].length();
2890 for (unsigned i = 0; i < nHits; i++) {
2891 TMLink & l = * _hits[2][i];
2892 links[l.wire()->superLayerId()].append(l);
2893 }
2894
2895 //...MC tracks...
2896 const unsigned nHep = TTrackHEP::list().length();
2897
2898 //...Loop over each super layer...
2899 for (unsigned i = 0; i < 11; i++) {
2900 unsigned superLayerId = i / 2;
2901 unsigned id = i % 2;
2902
2903 //...Counters...
2904 AList<TMLink> ** seeds =
2905 (AList<TMLink> **) malloc(nHep * sizeof(AList<TMLink> *));
2906 for (unsigned j = 0; j < nHep; j++)
2907 seeds[j] = NULL;
2908
2909 //...Hit loop...
2910 unsigned n = links[i].length();
2911 for (unsigned j = 0; j < n; j++) {
2912 TMLink & l = * links[i][j];
2913 const TMDCWireHitMC * mc = l.hit()->mc();
2914 if (! l.hit()->mc()) {
2915 std::cout << "TConformalFinder::findSegmentsPerfect !!! "
2916 << "no MC info. found ... aborted" << std::endl;
2917 return;
2918 }
2919 const unsigned id = mc->hep()->id();
2920
2921 if (! seeds[id])
2922 seeds[id] = new AList<TMLink>();
2923 seeds[id]->append(l);
2924 }
2925
2926 //...Create segments...
2927 AList<TSegment> segments;
2928 for (unsigned j = 0; j < nHep; j++) {
2929 if (! seeds[j]) continue;
2930
2931 const unsigned length = seeds[j]->length();
2932 if (length < _minNLinksForSegment) continue;
2933 if (length > _maxNLinksForSegment) continue;
2934 if (Width(* seeds[j]) > _maxWidthForSegment) continue;
2935
2936 TSegment & s = * new TSegment();
2937 s.append(* seeds[j]);
2938 segments.append(s);
2939 }
2940 _allSegments[id][superLayerId].append(segments);
2941 _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
2942
2943 //...Clear...
2944 for (unsigned j = 0; j < nHep; j++) {
2945 if (seeds[j])
2946 delete seeds[j];
2947 }
2948 free(seeds);
2949 }
2950
2951#ifdef TRKRECO_DEBUG_DETAIL
2952 std::cout << "... segment finding(perfect) finished" << std::endl;
2953#endif
2954}
2955
2956void
2957TConformalFinder::findSegmentsTsf(void) {
2958 //...Create lists of links for each super layer...
2959 AList<TMLink> links[11];
2960 unsigned n = _hits[2].length();
2961 for (unsigned i = 0; i < n; i++) {
2962 TMLink & l = * _hits[2][i];
2963 links[l.wire()->superLayerId()].append(l);
2964 }
2965
2966 //find Segments by tsf
2967 TMDCTsf * tsf[11];
2968
2969 //initial
2970 for (unsigned i = 0; i < 11; ++i)
2971 tsf[i] = new TMDCTsf(i);
2972
2973 for (unsigned i = 0; i < 11; i++) {
2974 unsigned superLayerId;
2975 unsigned id;
2976 switch(i){
2977 case 0: superLayerId = 0; id = 1; break;
2978 case 1: superLayerId = 1; id = 1; break;
2979 case 2: superLayerId = 0; id = 0; break;
2980 case 3: superLayerId = 1; id = 0; break;
2981 case 4: superLayerId = 2; id = 0; break;
2982 case 5: superLayerId = 2; id = 1; break;
2983 case 6: superLayerId = 3; id = 1; break;
2984 case 7: superLayerId = 4; id = 1; break;
2985 case 8: superLayerId = 5; id = 1; break;
2986 case 9: superLayerId = 3; id = 0; break;
2987 case 10: superLayerId = 4; id = 0; break;
2988 default: break;
2989 }
2990 //find segments
2991 AList<TSegment> tmp = tsf[i]->segments(links[i]);
2992
2993 //...Store them...
2994 _allSegments[id][superLayerId].append(tmp);
2995 _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
2996 delete tsf[i];
2997
2998 //...remove used in seg...
2999// for (int j = 0; j < _allUnused[id][superLayerId].length(); ++j){
3000// AList<TMLink> used = _allUnused[id][superLayerId][j]->links();
3001// _unused[id].remove(used);
3002// _unused[2].remove(used);
3003// }
3004 }
3005}
HepVector3D
HepGeom::Vector3D< double > HepVector3D
Definition Analysis/VertexFit/VertexFit-00-03-05/src/BField.cxx:13
n
const Int_t n
Definition DataBase/tau_mode.c:65
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
phi1
Double_t phi1
Definition DataBase/tau_mode.c:7
output
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in output
Definition FoamA.h:89
unit
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
Definition FoamA.h:90
s
XmlRpcServer s
Definition HelloServer.cpp:11
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition KarLud.h:35
n1
int n1
Definition SD0Tag.cxx:54
M_PI
#define M_PI
Definition TConstant.h:4
bitDisplay
void bitDisplay(unsigned val)
Definition TMDCUtil.cxx:85
ORIGIN
const HepPoint3D ORIGIN
Constants.
Definition TMDCUtil.cxx:47
WireHitConformalFinder
#define WireHitConformalFinder
Definition TMDCWireHit.h:51
WireHitUsed
#define WireHitUsed
Definition TMDCWireHit.h:47
NMajorLinks
unsigned NMajorLinks(const TSegment &a)
returns # of links in the major link.
Definition TSegment.cxx:1005
NUniqueLinks
unsigned NUniqueLinks(const TSegment &a)
checks property of segments.
Definition TSegment.cxx:978
OuterMostUniqueLink
TSegment * OuterMostUniqueLink(const TSegment &a)
returns a segment to the outer-most unique segment.
Definition TSegment.cxx:1084
UniqueLinks
AList< TSegment > UniqueLinks(const TSegment &a)
returns a list of unique segments in links.
Definition TSegment.cxx:991
MajorLinks
AList< TSegment > MajorLinks(const TSegment &a)
returns a list of segments in major links.
Definition TSegment.cxx:1026
Links
AList< TMLink > Links(const TSegment &s, const TTrack &t)
returns AList of TMLink used for a track.
Definition TSegment.cxx:963
TrackFastFinder
#define TrackFastFinder
Definition TTrack.h:24
TrackQuality2D
#define TrackQuality2D
Definition TTrack.h:47
TrackSlowFinder
#define TrackSlowFinder
Definition TTrack.h:25
TrackOldConformalFinder
#define TrackOldConformalFinder
Definition TTrack.h:23
TrackFitGlobal
#define TrackFitGlobal
Definition TTrack.h:52
NULL
#define NULL
Definition VertexExtrapolate.h:14
t
TTree * t
Definition binning.cxx:23
IMagneticFieldSvc::getReferField
virtual double getReferField()=0
TBuilder::buildRphi
TTrack * buildRphi(const AList< TMLink > &) const
builds a r/phi track from TMLinks or from Segments.
Definition TBuilder.cxx:107
TBuilder::salvage
void salvage(TTrack &t, AList< TMLink > &hits) const
salvages hits.
Definition TBuilder.cxx:211
TBuilder::buildStereoNew
TTrack * buildStereoNew(const TTrack &t, AList< TSegment > &goodSegments, AList< TSegment > &badSegments) const
Definition TBuilder.cxx:1535
TBuilder::buildStereo
TTrack * buildStereo(const TTrack &t, AList< TSegment > &) const
Definition TBuilder.cxx:658
TCircle
A class to represent a circle in tracking.
Definition TCircle.h:42
TConformalFinder0::conformalTransformationDriftCircle
static void conformalTransformationDriftCircle(const HepPoint3D &center, const AList< TMDCWireHit > &hits, AList< TMLink > &links)
transforms drift circle of hits into a conformal plane. transformed positions( x0,...
Definition TConformalFinder0.cxx:127
TConformalFinder::dump
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TConformalFinder.cxx:123
TConformalFinder::clear
void clear(void)
clear internal information.
Definition TConformalFinder.cxx:173
TConformalFinder::doit
int doit(const AList< TMDCWireHit > &axialHits, const AList< TMDCWireHit > &stereoHits, AList< TTrack > &tracks, AList< TTrack > &tracks2D)
finds tracks.
Definition TConformalFinder.cxx:806
TConformalFinder::version
std::string version(void) const
returns version.
Definition TConformalFinder.cxx:55
TConformalFinder::conformal2Wire
static const TMDCWire * conformal2Wire(const HepPoint3D &conformalPoint)
Definition TConformalFinder.cxx:1529
TConformalFinder::name
std::string name(void) const
returns name.
Definition TConformalFinder.h:262
TConformalFinder::~TConformalFinder
virtual ~TConformalFinder()
Destructor.
Definition TConformalFinder.cxx:119
TConformalFinder::TConformalFinder
TConformalFinder(unsigned fastFinder, unsigned slowFinder, unsigned perfectSegmentFinding, float maxSigma, float maxSigmaStereo, float salvageLevel, unsigned minNLinksForSegment, unsigned minNCoreLinks, unsigned minNSegments, unsigned salvageLoadWidth, unsigned stereoMode, unsigned stereoLoadWidth, float szSegmentDistance, float szLinkDistance, unsigned fittingFlag)
Constructor.
Definition TConformalFinder.cxx:59
TFinderBase::debugLevel
virtual int debugLevel(void) const
returns debug level.
Definition TFinderBase.h:90
TFinderBase::dump
virtual void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TFinderBase.cxx:23
THistogram
A class for a histogram used in tracking.
Definition THistogram.h:34
THistogram::segments
AList< TSegment > segments(void) const
returns an AList<TSegment0> using clusters() function.
Definition THistogram.cxx:259
THistogram::fillPhi
void fillPhi(const AList< TMLink > &links)
fills with hits.
Definition THistogram.cxx:111
TMDCTsf
A class to represent a Track Finder Segment(TSF).
Definition TMDCTsf.h:35
TMDCTsf::segments
AList< TSegment > segments(const AList< TMLink > &links)
finds segments.
Definition TMDCTsf.cxx:199
TMDCWireHitMC
A class to represent a MC wire hit in MDC.
Definition TMDCWireHitMC.h:46
TMDCWireHitMC::hep
const TTrackHEP *const hep(void) const
returns a pointer to a GEN_HEPEVT.
Definition TMDCWireHitMC.h:159
TMDCWireHit::drift
float drift(unsigned) const
returns drift distance.
Definition TMDCWireHit.h:236
TMDCWireHit::mc
const TMDCWireHitMC *const mc(void) const
returns a pointer to TMDCWireHitMC.
Definition TMDCWireHit.h:298
TMDCWire
A class to represent a wire in MDC.
Definition TMDCWire.h:55
TMDCWire::superLayerId
unsigned superLayerId(void) const
returns super layer id.
Definition TMDCWire.h:225
TMDCWire::name
std::string name(void) const
returns name.
Definition TMDCWire.h:412
TMDC
Definition TMDC.h:61
TMDC::wire
const TMDCWire *const wire(unsigned wireId) const
returns a pointer to a wire. 0 will be returned if 'wireId' is invalid.
Definition TMDC.h:277
TMDC::getTMDC
static TMDC * getTMDC(void)
Definition TMDC.cxx:95
TMLine
A class to represent a track in tracking.
Definition TMLine.h:40
TPoint2D
A class to represent a point in 2D.
Definition TPoint2D.h:37
TPoint2D::mag2
double mag2(void) const
Definition TPoint2D.h:118
TPoint2D::cross
double cross(const TPoint2D &) const
Definition TPoint2D.h:139
TPoint2D::phi
double phi(void) const
Definition TPoint2D.h:124
TPoint2D::y
double y(void) const
Definition TPoint2D.h:94
TPoint2D::x
double x(void) const
Definition TPoint2D.h:88
TSegment
A class to relate TMDCWireHit and TTrack objects.
Definition TSegment.h:43
TSegment::direction
const HepVector3D & direction(void) const
returns direction.
Definition TSegment.h:284
TSegment::dump
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TSegment.cxx:116
TSegment::outerLinks
AList< TSegment > & outerLinks(void)
Definition TSegment.h:378
TSegment::position
const HepPoint3D & position(void) const
returns position.
Definition TSegment.h:277
TSegment::distance
double distance(const TSegment &) const
calculates distance between two clusters. Smaller value indicates closer.
Definition TSegment.cxx:237
TSegment::innerLinks
AList< TSegment > & innerLinks(void)
Definition TSegment.h:366
TSegment::tracks
AList< TTrack > & tracks(void)
Definition TSegment.h:360
TTrackBase
A virtual class for a track class in tracking.
Definition TTrackBase.h:46
TTrackBase::fit
virtual int fit(void)
fits itself by a default fitter. Error was happened if return value is not zero.
Definition TTrackBase.cxx:357
TTrackBase::links
const AList< TMLink > & links(unsigned mask=0) const
returns a list of masked TMLinks assigned to this track. 'mask' will be applied if mask is not 0.
Definition TTrackBase.cxx:297
TTrackBase::remove
void remove(TMLink &a)
removes a TMLink.
Definition TTrackBase.h:204
TTrackBase::nCores
unsigned nCores(unsigned mask=0) const
returns # of masked TMLinks for fit. 'mask' will be applied if mask is not 0.
Definition TTrackBase.cxx:325
TTrackBase::cores
const AList< TMLink > & cores(unsigned mask=0) const
returns a list of masked TMLinks for fit. 'mask' will be applied if mask is not 0.
Definition TTrackBase.cxx:317
TTrackHEP::id
unsigned id(void) const
returns an id started from 0.
Definition TTrackHEP.h:122
TTrackHEP::list
static const AList< TTrackHEP > & list(void)
returns a list of TTrackHEP's.
Definition TTrackHEP.cxx:72
TTrackManager::maskBadHits
static void maskBadHits(const AList< TTrack > &, float maxSigma2)
masks hits with large chisq as associated hits. Pull in TMLink is used.
Definition TTrackManager.cxx:2487
TTrackManager::goodTrack
static bool goodTrack(const TTrack &, bool track2D=false)
checks goodness of a track.
Definition TTrackManager.cxx:3015
TTrack
A class to represent a track in tracking.
Definition TTrack.h:129
TTrack::segments
AList< TSegment > & segments(void)
returns AList<TSegment>.
Definition TTrack.h:583
TTrack::name
const std::string & name(void) const
returns/sets name.
Definition TTrack.h:516
TTrack::center
TPoint2D center(void) const
returns position of helix center.
Definition TTrack.h:595
TTrack::radius
double radius(void) const
returns signed radius.
Definition TTrack.h:577
TTrack::approach
int approach(TMLink &) const
calculates the closest approach to a wire in real space. Results are stored in TMLink....
Definition TTrack.cxx:296
TTrack::charge
double charge(void) const
returns charge.
Definition TTrack.h:504
y
double y[1000]
Definition draw_charge_space_coarse_II.cxx:8
Alignment::PI2
const double PI2
Definition Alignment.h:41
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:27
b
const double b
Definition slope.cxx:9
NS
#define NS(x)
Definition xmltok.c:1503