BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
TBuilder.cxx
Go to the documentation of this file.
1//-----------------------------------------------------------------------------
2// $Id: TBuilder.cxx,v 1.28 2022/01/28 22:04:42 maqm Exp $
3//-----------------------------------------------------------------------------
4// Filename : TBuilder.cc
5// Section : Tracking
6// Owner : Yoshi Iwasaki
7// Email : [email protected]
8//-----------------------------------------------------------------------------
9// Description : A class to build a track.
10// See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
11//-----------------------------------------------------------------------------
12
13#include "CLHEP/String/Strings.h"
14#include "TrkReco/TBuilder.h"
15#include "TrkReco/TMDC.h"
16#include "TrkReco/TMDCWire.h"
17#include "TrkReco/TMDCWireHit.h"
18#include "TrkReco/TMLink.h"
19#include "TrkReco/TCircle.h"
20#include "TrkReco/TLine0.h"
21#include "TrkReco/TMLine.h"
22#include "TrkReco/TTrack.h"
23#include "TrkReco/TSegment.h"
24#include "TrkReco/TPoint2D.h"
25#include "TrkReco/TLine2D.h"
26#include "TrkReco/TRobustLineFitter.h"
27#ifdef TRKRECO_DEBUG
28#include "TrkReco/TMDCWireHitMC.h"
29#include "TrkReco/TTrackHEP.h"
30#endif
31#ifdef TRKRECO_WINDOW
32#include "TrkReco/TWindow.h"
33TWindow sz("sz");
34#endif
35
36TBuilder::TBuilder(const std::string & a,
37 float maxSigma,
38 float maxSigmaStereo,
39 float salvageLevel,
40 float szSegmentDistance,
41 float szLinkDistance,
42 unsigned fittingFlag)
43: _name(a),
44 _fitter("TBuilder Fitter"),
45 _maxSigma(maxSigma),
46 _maxSigmaStereo(maxSigmaStereo),
47 _salvageLevel(sqrt(salvageLevel)),
48// _minNCores(3), //20060307
49 _minNCores(2), //liuqg20070312
50 _szSegmentDistance(szSegmentDistance),
51 _szLinkDistance(szLinkDistance) {
52 if (fittingFlag & 1) _fitter.sag(true);
53 if (fittingFlag & 2) _fitter.propagation(true);
54 if (fittingFlag & 4) _fitter.tof(true);
55 if (fittingFlag & 8) _fitter.freeT0(true);
56}
57
60
61void
62TBuilder::dump(const std::string & msg, const std::string & pre) const {
63}
64
65TTrack *
66TBuilder::buildRphi(const AList<TSegment> & list) const {
67
68#ifdef TRKRECO_DEBUG
69 std::cout << "... building rphi by segments : # of segments = ";
70 std::cout << list.length() << std::endl;
71 for (unsigned i = 0; i < list.length(); i++)
72 list[i]->dump("hits sort flag", " ");
73#endif
74 //yuany
75 // for (unsigned i = 0; i < list.length(); i++)
76 // list[i]->dump("", " ");
77
78 //...Pick up links...
79 AList<TMLink> links = Links(list);
80 //yuany
81// for(unsigned i=0;i<links.length();i++){
82// cout<<"wire name "<<links[i]->wire()->name()<<endl;
83// }
84
85 //...Main funtion...
86 TTrack * t = buildRphi(links);
87
88 //...Check used segments...
89 if (t) {
90 const AList<TMLink> & usedLinks = t->links();
91 unsigned n = list.length();
92 for (unsigned i = 0; i < n; i++) {
93 TSegment & segment = * list[i];
94 AList<TMLink> used = Links(segment, * t);
95 if (used.length()) {
96 t->segments().append(segment);
97 segment.tracks().append(t);
98 }
99 }
100 }
101
102 return t;
103
104}
105
106TTrack *
107TBuilder::buildRphi(const AList<TMLink> & list) const {
108#ifdef TRKRECO_DEBUG
109 std::cout << "... building rphi by links : # of links = ";
110 std::cout << list.length() << std::endl;
111#endif
112
113 //...Classify TMLink's...
114 AList<TMLink> cores;
115 AList<TMLink> nonCores;
116 SeparateCores(list, cores, nonCores);
117
118#ifdef TRKRECO_DEBUG
119 cout<<" ... cores..."<<endl;
120 for(int ii = 0; ii < cores.length(); ++ii) {
121 cout<<"layer: "<<cores[ii]->wire()->layerId()
122 <<" local: "<<cores[ii]->wire()->localId()<<endl;
123 }
124 cout<<" ...noncores..."<<endl;
125 for(int ii = 0; ii < nonCores.length(); ++ii) {
126 cout<<"layer: "<<nonCores[ii]->wire()->layerId()
127 <<" local: "<<nonCores[ii]->wire()->localId()<<endl;
128 }
129#endif
130 //...Check # of links...
131 unsigned nCores = cores.length();
132 if (nCores < _minNCores) {
133#ifdef TRKRECO_DEBUG
134 std::cout << "... building rphi failure : # of cores(=" << nCores;
135 std::cout << ") is less then " << _minNCores << std::endl;
136#endif
137
138 return NULL;
139 }
140
141 //...Make a circle...
142#ifdef TRKRECO_DEBUG
143 std::cout <<"links in list = "<<list.length()<<std::endl;
144 std::cout << "... making a circle : # cores =" << cores.length() << std::endl;
145#endif
146 TCircle c(cores);
147 int err = c.fit();
148 if (err < 0) {
149#ifdef TRKRECO_DEBUG
150 std::cout << "... building rphi failure : circle fit error = ";
151 std::cout << err << std::endl;
152#endif
153 return NULL;
154 }
155// cout<<"radius "<<c.radius()<<" charge "<<c.charge()<<" center "<<c.center()<<endl;
156 //...Make a track...
157 TTrack * t = new TTrack(c);
158// t->dump("hits sort flag", " ");
159 AList<TMLink> bad;
160 err = _fitter.fit(* t);
161 //cout<<"fit 1 done err: "<<err<<endl;
162 // t->dump("hits sort flag", " ");
163 if (err < 0) goto discard;
164 t->refine(bad, _maxSigma * 100.);
165 //cout<<"refine 1 done"<<endl;
166 //t->dump("hits sort flag", " ");
167 err = _fitter.fit(* t);
168 //cout<<"fit 2 done err: "<<err<<endl;
169 //t->dump("hits sort flag", " ");
170 t->refine(bad, _maxSigma * 10.);
171 //cout<<"refine 2 done"<<endl;
172 //t->dump("hits sort flag", " ");
173 err = _fitter.fit(* t);
174 //cout<<"fit 3 done err: "<<err<<endl;
175 //t->dump("hits sort flag", " ");
176 t->refine(bad, _maxSigma);
177 //cout<<"refine 3 done"<<endl;
178 //t->dump("hits sort flag", " ");
179 err = _fitter.fit(* t);
180 //cout<<"fit 4 done err: "<<err<<endl;
181
182 if (err < 0) goto discard;
183#ifdef TRKRECO_DEBUG_DETAIL
184 c.dump("detail", " ccl> ");
185 t->dump("detail", " 1st> ");
186#endif
187
188 //...Try to append non-core hits...
189#ifdef TRKRECO_DEBUG
190 std::cout << "... appending non-core hits : # = " << nonCores.length();
191 std::cout << std::endl;
192#endif
193 t->appendByApproach(nonCores, _salvageLevel);
194#ifdef TRKRECO_DEBUG
195 t->dump("hits sort flag", " ");
196#endif
197// t->dump("hits sort flag", " ");
198 return t;
199
200 //...Something happened...
201discard:
202#ifdef TRKRECO_DEBUG
203 std::cout << "... building rphi failure : helix fit error = ";
204 std::cout << err << std::endl;
205#endif
206 delete t;
207 return NULL;
208}
209
210void
211TBuilder::salvage(TTrack & t, AList<TMLink> & hits) const {
212
213#ifdef TRKRECO_DEBUG
214 std::cout << "... salvaging(TBuilder) : # of given hits=" << hits.length();
215 std::cout << ", salvage level=" << _salvageLevel << std::endl;
216 Dump(hits, "hits sort flag");
217#endif
218
219 unsigned nHits = hits.length();
220 if (nHits == 0) return;
221
222 //...Try to append this hit...
223 t.appendByApproach(hits, _salvageLevel);
224 _fitter.fit(t);
225}
226
227void
228TBuilder::salvage(TTrack & t, AList<TSegment> & list) const {
229
230#ifdef TRKRECO_DEBUG
231 std::cout << "... salvaging by segments : # of segments = ";
232 std::cout << list.length() << std::endl;
233 for (unsigned i = 0; i < list.length(); i++)
234 list[i]->dump("hits sort flag", " ");
235#endif
236
237 //...Pick up links...
238 AList<TMLink> links;
239 unsigned n = list.length();
240 for (unsigned i = 0; i < n; i++)
241 links.append(list[i]->links());
242
243 salvage(t, links);
244}
245
246TMLine *
247TBuilder::initialLine1(const TTrack & t,
248 const AList<TSegment> & list,
249 const AList<TMLink> & lList) const {
250
251 //...Check input...
252 if (list.length() < 2) return NULL;
253
254 //...Check super layer pattern...
255 AList<TSegment> sl[5];
256 AList<TMLink> tl[5];
257 unsigned n = list.length();
258 for (unsigned i = 0; i < n; i ++) {
259 unsigned j = list[i]->superLayerId() / 2;
260 sl[j].append(list[i]);
261 tl[j].append(lList[i]);
262 }
263#ifdef TRKRECO_DEBUG
264 std::cout << " ... initialLine1 : super layer ptn = ";
265 for (unsigned i = 0; i < 5; i++) std::cout << sl[i].length();
266 std::cout << std::endl;
267#endif
268
269 //...Count single segment layer...
270 unsigned nSingle = 0;
271 for (unsigned i = 0; i < 5; i++)
272 if (sl[i].length() == 1) ++nSingle;
273#ifdef TRKRECO_DEBUG
274 std::cout << " ... # of single segment layer = " << nSingle << std::endl;
275#endif
276 if (nSingle < 2) return NULL;
277
278 //...Happy case...
279 AList<TSegment> bestCombination;
280 AList<TMLink> forLine;
281 for (unsigned i = 0; i < 5; i++) {
282 if (sl[i].length() != 1) continue;
283 bestCombination.append(sl[i]);
284 forLine.append(tl[i]);
285 }
286 TMLine & line = * new TMLine(forLine);
287 int err = line.fit();
288 if (err) {
289 delete & line;
290 return NULL;
291 }
292 return & line;
293}
294
295TMLine *
296TBuilder::initialLine2(const TTrack & t,
297 const AList<TMLink> & lList) const {
298#ifdef TRKRECO_DEBUG
299 std::cout << " ... initlialLine2 : # of links = " << lList.length();
300 std::cout << std::endl;
301#endif
302
303 //...Static object...
304 static TRobustLineFitter fitter("can you work?");
305 TMLine * line = new TMLine(lList);
306 int err = fitter.fit(* line);
307 if (err) {
308 delete line;
309 return NULL;
310 }
311 return line;
312}
313
314TMLine *
315TBuilder::initialLine(const TTrack & t, AList<TSegment> & list) const {
316
317 //...Select good segments in acceptance...
318 AList<TSegment> bad;
319 AList<TMLink> lList;
320 unsigned n = list.length();
321 for (unsigned i = 0; i < n; i++) {
322 TMLink & l = * new TMLink();
323 int err = t.szPosition(* list[i], l);
324 if (err) {
325 bad.append(list[i]);
326 delete & l;
327 continue;
328 }
329
330 const HepPoint3D & p = l.position();
331 if (p.y() > l.wire()->forwardPosition().z() ||
332 p.y() < l.wire()->backwardPosition().z()) {
333 bad.append(list[i]);
334 delete & l;
335 continue;
336 }
337 lList.append(l);
338 }
339 list.remove(bad);
340 n = list.length();
341
342 //...Cal. expected # of super layers...
343 unsigned nAMaxSL = 0;
344 unsigned nAMinSL = 10;
345 const AList<TSegment> & segments = t.segments();
346 unsigned nA = segments.length();
347 for (unsigned i = 0; i < nA; i++) {
348 unsigned sl = segments[i]->superLayerId();
349 if (sl > nAMaxSL) nAMaxSL = sl;
350 if (sl < nAMinSL) nAMinSL = sl;
351 }
352 unsigned nExpected = (nAMaxSL - nAMinSL) / 2;
353#ifdef TRKRECO_DEBUG
354 std::cout << " ... initialLine : axial super layer usage = " << nAMinSL;
355 std::cout << " ~ " << nAMaxSL << std::endl;
356 std::cout << " : expected stereo super layers = ";
357 std::cout << nExpected << std::endl;
358#endif
359
360 TMLine * line = NULL;
361 AList<TSegment> list0 = list;
362 AList<TMLink> lList0 = lList;
363 while (1) {
364 bool last = (NSuperLayers(lList) <= nExpected);
365
366 TMLine * line = initialLine2(t, lList);
367 if (line) {
368 AList<TSegment> tmp = selectStereoSegment(* line, list0, lList0);
369#ifdef TRKRECO_DEBUG
370 std::cout << " ... initialLine2 : # of selected segments = ";
371 std::cout << tmp.length() << std::endl;
372#endif
373 if ((tmp.length() >= nExpected) ||
374 (tmp.length() >= 4) ||
375 last) {
376 list = tmp;
377 return line;
378 }
379 delete line;
380 }
381
382 line = initialLine1(t, list, lList);
383 if (line) {
384 AList<TSegment> tmp = selectStereoSegment(* line, list0, lList0);
385 bool ok = (tmp.length() >= nExpected);
386#ifdef TRKRECO_DEBUG
387 std::cout << " ... initialLine2 : # of selected segments = ";
388 std::cout << tmp.length() << " : ok, last = ";
389 std::cout << ok << ", " << last << std::endl;
390#endif
391 if (ok) {
392 list = tmp;
393 return line;
394 }
395 removeFarSegment(* line, list, lList);
396 delete line;
397 }
398 else {
399 return NULL;
400 }
401 }
402}
403
404TMLine *
405TBuilder::initialLineOld(const TTrack & t, AList<TSegment> & list) const {
406
407 //...Select good segments in acceptance...
408 AList<TSegment> bad;
409 AList<TMLink> lList;
410 unsigned n = list.length();
411 for (unsigned i = 0; i < n; i++) {
412 TMLink & l = * new TMLink();
413 int err = t.szPosition(* list[i], l);
414 if (err) {
415 bad.append(list[i]);
416 delete & l;
417 continue;
418 }
419
420 const HepPoint3D & p = l.position();
421 if (p.y() > l.wire()->forwardPosition().z() ||
422 p.y() < l.wire()->backwardPosition().z()) {
423 bad.append(list[i]);
424 delete & l;
425 continue;
426 }
427 lList.append(l);
428 }
429 list.remove(bad);
430 n = list.length();
431
432 //...Cal. expected # of super layers...
433 unsigned nAMaxSL = 0;
434 unsigned nAMinSL = 10;
435 const AList<TSegment> & segments = t.segments();
436 unsigned nA = segments.length();
437 for (unsigned i = 0; i < nA; i++) {
438 unsigned sl = segments[i]->superLayerId();
439 if (sl > nAMaxSL) nAMaxSL = sl;
440 if (sl < nAMinSL) nAMinSL = sl;
441 }
442 unsigned nExpected = (nAMaxSL - nAMinSL) / 2;
443#ifdef TRKRECO_DEBUG
444 std::cout << " ... initialLine : axial super layer usage = " << nAMinSL;
445 std::cout << " ~ " << nAMaxSL << std::endl;
446 std::cout << " : expected stereo super layers = ";
447 std::cout << nExpected << std::endl;
448#endif
449
450 TMLine * line = NULL;
451 AList<TSegment> list0 = list;
452 AList<TMLink> lList0 = lList;
453 while (1) {
454 TMLine * line = initialLine1(t, list, lList);
455 if (line) {
456 AList<TSegment> tmp = selectStereoSegment(* line, list0, lList0);
457#ifdef TRKRECO_DEBUG
458 std::cout << " ... initialLine1 : # of selected segments = ";
459 std::cout << tmp.length() << std::endl;
460#endif
461 if (tmp.length() >= nExpected) {
462 lList0.remove(line->links());
463 HepAListDeleteAll(lList0);
464 list = tmp;
465 return line;
466 }
467 delete line;
468 }
469
470 line = initialLine2(t, lList);
471 if (line) {
472 AList<TSegment> tmp = selectStereoSegment(* line, list0, lList0);
473 bool ok = (tmp.length() >= nExpected);
474 bool last = (NSuperLayers(lList) <= nExpected);
475#ifdef TRKRECO_DEBUG
476 std::cout << " ... initialLine2 : # of selected segments = ";
477 std::cout << tmp.length() << " : ok, last = ";
478 std::cout << ok << ", " << last << std::endl;
479#endif
480 if (ok || last) {
481 lList0.remove(line->links());
482 HepAListDeleteAll(lList0);
483 list = tmp;
484 return line;
485 }
486 removeFarSegment(* line, list, lList);
487 delete line;
488 }
489 else {
490 HepAListDeleteAll(lList0);
491 return NULL;
492 }
493 }
494}
495
496AList<TSegment>
497TBuilder::selectStereoSegment(const TMLine & line,
498 const AList<TSegment> & list,
499 const AList<TMLink> & szList) const {
500 AList<TSegment> outList;
501 unsigned n = list.length();
502 for (unsigned i = 0; i < n; i++) {
503 float distance = line.distance(* szList[i]);
504 if (distance < _szSegmentDistance) outList.append(list[i]);
505 }
506 return outList;
507}
508
509void
510TBuilder::removeFarSegment(const TMLine & line,
511 AList<TSegment> & list,
512 AList<TMLink> & szList) const {
513 unsigned n = list.length();
514 float maxDistance = 0.;
515 unsigned maxId = 0;
516 for (unsigned i = 0; i < n; i++) {
517 float distance = line.distance(* szList[i]);
518 if (distance > maxDistance) maxId = i;
519 }
520 list.remove(maxId);
521 szList.remove(maxId);
522}
523
524TTrack *
525TBuilder::buildStereo(TTrack & track,
526 TMLine & line,
527 const AList<TMLink> & list) const {
528#ifdef TRKRECO_DEBUG
529 std::cout << "... building stereo by links : # of links = ";
530 std::cout << list.length() << std::endl;
531#endif
532
533 //...Classify TMLink's...
534 AList<TMLink> cores;
535 AList<TMLink> nonCores;
536 SeparateCores(list, cores, nonCores);
537 //...Check # of links...
538 unsigned nCores = cores.length();
539 if (nCores < _minNCores) {
540#ifdef TRKRECO_DEBUG
541 std::cout << "... stereo building failure : # of cores(=" << nCores;
542 std::cout << ") is less then " << _minNCores << std::endl;
543#endif
544 return NULL;
545 }
546#ifdef TRKRECO_WINDOW
547 sz.appendSz(track, cores, leda_brown);
548#endif
549
550 //...Cal. left and right position...
551 AList<TMLink> forNewLine;
552 for (unsigned i = 0; i < nCores; i++) {
553 TMLink & t = * cores[i];
554 for (unsigned i = 0; i < 2; i++) {
555 TMLink & tt = * new TMLink(t);
556 tt.leftRight(i);
557 int err = track.szPosition(tt);
558 if (err) {
559 delete & tt;
560 }
561 else {
562 if (line.distance(tt) < _szLinkDistance) {
563 tt.link(& t);
564 forNewLine.append(tt);
565 }
566 else {
567 delete & tt;
568 }
569 }
570 }
571 }
572
573 //...Create new line...
574#ifdef TRKRECO_DEBUG
575 std::cout << " ... creating a new line" << std::endl;
576#endif
577 unsigned nNewLine = forNewLine.length();
578 TMLine newLine(forNewLine);
579 static TRobustLineFitter fitter("abc");
580 int err = fitter.fit(newLine);
581 // int err = newLine.fit();
582 if (err < 0) {
583 HepAListDeleteAll(forNewLine);
584#ifdef TRKRECO_DEBUG_DETAIL
585 std::cout << " ... 2nd linear fit failure. nLinks(";
586 std::cout << forNewLine.length() << ")" << std::endl;
587#endif
588 return NULL;
589 }
590
591#ifdef TRKRECO_WINDOW
592 sz.append(newLine, leda_green);
593 sz.wait();
594#endif
595
596#ifdef TRKRECO_DEBUG_DETAIL
597 Dump(forNewLine, "sort hits stereo", " ");
598#endif
599
600 //...Decide left/right...
601 AList<TMLink> toRemove;
602 TMLink * last = NULL;
603 for (unsigned i = 0; i < nNewLine; i++) {
604 if (last == NULL) last = forNewLine[i]->link();
605 else {
606 if (last == forNewLine[i]->link()) {
607 if (newLine.distance(* forNewLine[i - 1]) >
608 newLine.distance(* forNewLine[i]))
609 toRemove.append(forNewLine[i - 1]);
610 else
611 toRemove.append(forNewLine[i]);
612 last = NULL;
613 }
614 else {
615 last = forNewLine[i]->link();
616 }
617 }
618 }
619 forNewLine.remove(toRemove);
620 nNewLine = forNewLine.length();
621
622#ifdef TRKRECO_DEBUG_DETAIL
623 Dump(toRemove, "sort hits stereo", " x ");
624 Dump(forNewLine, "sort hits stereo", " ");
625#endif
626
627 //...3D fit...
628 for (unsigned i = 0; i < nNewLine; i++)
629 track.append(* forNewLine[i]->link());
630 Vector a(5);
631 a = track.helix().a();
632 a[3] = newLine.b();
633 a[4] = track.charge() * newLine.a();
634 track._helix->a(a);
635
636 //...Refine...
637 AList<TMLink> bad;
638 err = _fitter.fit(track);
639 track.refine(bad, _maxSigmaStereo * 100.);
640 err = _fitter.fit(track);
641 track.refine(bad, _maxSigmaStereo * 10.);
642 err = _fitter.fit(track);
643 track.refine(bad, _maxSigmaStereo);
644 err = _fitter.fit(track);
645
646#ifdef TRKRECO_WINDOW
647 sz.text("stereo finished");
648 sz.oneShot(track, leda_blue);
649#endif
650
651 //...Termination...
652 HepAListDeleteAll(toRemove);
653 HepAListDeleteAll(forNewLine);
654 return & track;
655}
656
657TTrack *
658TBuilder::buildStereo(const TTrack & t0, AList<TSegment> & segments) const {
659
660 TTrack & t = * new TTrack(t0);
661
662#ifdef TRKRECO_DEBUG
663 std::cout << "... building stereo by segments : # of segments = ";
664 std::cout << segments.length() << std::endl;
665 for (unsigned i = 0; i < segments.length(); i++)
666 segments[i]->dump("hits sort flag", " ");
667#endif
668#ifdef TRKRECO_WINDOW
669 sz.clear();
670 sz.skip(false);
671 sz.mode(2);
672 sz.appendSz(t, segments, leda_black);
673 AList<TSegment> tmps = segments;
674 std::string s = "# of segments = " + itostring(int(segments.length()));
675#endif
676
677 //...Find initial line...
678 TMLine * line = initialLineOld(t, segments);
679 if (! line) {
680#ifdef TRKRECO_DEBUG
681 std::cout << "... building stereo failure : no initial line found" << std::endl;
682#endif
683#ifdef TRKRECO_WINDOW
684 s = "no initial line found : " + s;
685 sz.text(s);
686 sz.wait();
687#endif
688 return NULL;
689 }
690#ifdef TRKRECO_WINDOW
691 sz.append(* line, leda_red);
692 sz.text(s);
693 sz.wait();
694#endif
695
696 //...Pick up links...
697 AList<TMLink> links;
698 unsigned n = segments.length();
699 for (unsigned i = 0; i < n; i++)
700 links.append(segments[i]->links());
701
702 //...Main funtion...
703 TTrack * ts = buildStereo(t, * line, links);
704
705 //...Check used segments...
706 if (ts) {
707 AList<TMLink> usedLinks = StereoHits(t.links());
708 for (unsigned i = 0; i < n; i++) {
709 TSegment & segment = * segments[i];
710 AList<TMLink> used = Links(segment, t);
711 if (used.length()) {
712 t.segments().append(segment);
713 segment.tracks().append(t);
714 }
715 }
716 }
717
718 HepAListDeleteAll((AList<TMLink> &) line->links());
719 delete line;
720 return ts;
721}
722
723AList<TMLine>
724TBuilder::searchInitialLines(unsigned nSuperLayerMax) const {
725#ifdef TRKRECO_DEBUG
726 std::cout << " ... searchInitialLines : # of segments in sl : ";
727 for (unsigned i = 0; i < 5; i++)
728 std::cout << _links[i].length() << ",";
729 std::cout << std::endl
730 << " : max # of super layers="
731 << nSuperLayerMax << std::endl;
732#endif
733
734 AList<TMLine> lines;
735 if (nSuperLayerMax > 5)
736 lines.append(searchLines6());
737 else if (nSuperLayerMax > 4)
738 lines.append(searchLines5());
739 else if (nSuperLayerMax > 3)
740 lines.append(searchLines4());
741 else if (nSuperLayerMax > 2)
742 lines.append(searchLines3());
743 else if (nSuperLayerMax > 1)
744 lines.append(searchLines2());
745// else
746// lines.append(searchLines1());
747
748 lines.sort(SortByB);
749 return lines;
750}
751
752AList<TMLine>
753TBuilder::searchLines6(void) const {
754 unsigned n[6];
755 for (unsigned i = 0; i < 6; i++)
756 n[i] = _links[i].length();
757
758 AList<TMLine> lines;
759 for (unsigned i0 = 0; i0 < n[0]; i0++) {
760 for (unsigned i1 = 0; i1 < n[1]; i1++) {
761 for (unsigned i2 = 0; i2 < n[2]; i2++) {
762 for (unsigned i3 = 0; i3 < n[3]; i3++) {
763 for (unsigned i4 = 0; i4 < n[4]; i4++) {
764 for (unsigned i5 = 0; i5 < n[5]; i5++) {
765
766 AList<TMLink> forLine;
767 forLine.append(_links[0][i0]);
768 forLine.append(_links[1][i1]);
769 forLine.append(_links[2][i2]);
770 forLine.append(_links[3][i3]);
771 forLine.append(_links[4][i4]);
772 forLine.append(_links[5][i5]);
773
774 TMLine & line = * new TMLine(forLine);
775 int err = line.fit();
776 if (err) {
777 delete & line;
778 continue;
779 }
780
781 lines.append(line);
782 }
783 }
784 }
785 }
786 }
787 }
788
789#ifdef TRKRECO_DEBUG
790 std::cout << " ... searchLines5 : # of lines found = "
791 << lines.length() << std::endl;
792#endif
793
794 return lines;
795}
796
797AList<TMLine>
798TBuilder::searchLines5(void) const {
799 AList<TMLine> lines;
800
801 const AList<TMLink> * l[5];
802 for (unsigned skip = 0; skip < 6; skip++) {
803 if (skip == 0) {
804 l[0] = & _links[1];
805 l[1] = & _links[2];
806 l[2] = & _links[3];
807 l[3] = & _links[4];
808 l[4] = & _links[5];
809 }
810 else if (skip == 1) {
811 l[0] = & _links[0];
812 }
813 else if (skip == 2) {
814 l[1] = & _links[1];
815 }
816 else if (skip == 3) {
817 l[2] = & _links[2];
818 }
819 else if (skip == 4) {
820 l[3] = & _links[3];
821 }
822 else if (skip == 5) {
823 l[4] = & _links[4];
824 }
825
826 unsigned n[5];
827 for (unsigned i = 0; i < 5; i++)
828 n[i] = l[i]->length();
829
830 for (unsigned i0 = 0; i0 < n[0]; i0++) {
831 for (unsigned i1 = 0; i1 < n[1]; i1++) {
832 for (unsigned i2 = 0; i2 < n[2]; i2++) {
833 for (unsigned i3 = 0; i3 < n[3]; i3++) {
834 for (unsigned i4 = 0; i4 < n[4]; i4++) {
835
836 AList<TMLink> forLine;
837 forLine.append((* l[0])[i0]);
838 forLine.append((* l[1])[i1]);
839 forLine.append((* l[2])[i2]);
840 forLine.append((* l[3])[i3]);
841 forLine.append((* l[4])[i4]);
842
843 TMLine & line = * new TMLine(forLine);
844 int err = line.fit();
845 if (err) {
846 delete & line;
847 continue;
848 }
849
850 lines.append(line);
851 }
852 }
853 }
854 }
855 }
856 }
857
858#ifdef TRKRECO_DEBUG
859 std::cout << " ... searchLines5 : # of lines found = "
860 << lines.length() << std::endl;
861#endif
862
863 return lines;
864}
865
866AList<TMLine>
867TBuilder::searchLines4(void) const {
868 AList<TMLine> lines;
869
870 const AList<TMLink> * l[4];
871 for (unsigned skip = 0; skip < 15; skip++) {
872 if (skip == 0) {
873 l[0] = & _links[0];
874 l[1] = & _links[1];
875 l[2] = & _links[2];
876 l[3] = & _links[3];
877 }
878 else if (skip == 1) {
879 l[3] = & _links[4];
880 }
881 else if (skip == 2) {
882 l[3] = & _links[5];
883 }
884 else if (skip == 3) {
885 l[2] = & _links[3];
886 l[3] = & _links[4];
887 }
888 else if (skip == 4) {
889 l[3] = & _links[5];
890 }
891 else if (skip == 5) {
892 l[2] = & _links[4];
893 }
894 else if (skip == 6) {
895 l[1] = & _links[2];
896 l[2] = & _links[3];
897 l[3] = & _links[4];
898 }
899 else if (skip == 7) {
900 l[3] = & _links[5];
901 }
902 else if (skip == 8) {
903 l[2] = & _links[4];
904 }
905 else if (skip == 9) {
906 l[1] = & _links[3];
907 }
908 else if (skip == 10) {
909 l[0] = & _links[1];
910 l[1] = & _links[2];
911 l[2] = & _links[3];
912 l[3] = & _links[4];
913 }
914 else if (skip == 11) {
915 l[3] = & _links[5];
916 }
917 else if (skip == 12) {
918 l[2] = & _links[4];
919 }
920 else if (skip == 13) {
921 l[1] = & _links[3];
922 }
923 else if (skip == 14) {
924 l[0] = & _links[2];
925 }
926
927 unsigned n[4];
928 for (unsigned i = 0; i < 4; i++)
929 n[i] = l[i]->length();
930
931 for (unsigned i0 = 0; i0 < n[0]; i0++) {
932 for (unsigned i1 = 0; i1 < n[1]; i1++) {
933 for (unsigned i2 = 0; i2 < n[2]; i2++) {
934 for (unsigned i3 = 0; i3 < n[3]; i3++) {
935 AList<TMLink> forLine;
936 forLine.append((* l[0])[i0]);
937 forLine.append((* l[1])[i1]);
938 forLine.append((* l[2])[i2]);
939 forLine.append((* l[3])[i3]);
940
941 TMLine & line = * new TMLine(forLine);
942 int err = line.fit();
943 if (err) {
944 delete & line;
945 continue;
946 }
947
948 lines.append(line);
949 }
950 }
951 }
952 }
953 }
954
955#ifdef TRKRECO_DEBUG
956 std::cout << " ... searchLines4 : # of lines found = "
957 << lines.length() << std::endl;
958#endif
959
960 return lines;
961}
962
963AList<TMLine>
964TBuilder::searchLines3(void) const {
965 AList<TMLine> lines;
966
967 const AList<TMLink> * l[3];
968 for (unsigned skip = 0; skip < 20; skip++) {
969 if (skip == 0) {
970 l[0] = & _links[0];
971 l[1] = & _links[1];
972 l[2] = & _links[2];
973 }
974 else if (skip == 1) {
975 l[2] = & _links[3];
976 }
977 else if (skip == 2) {
978 l[2] = & _links[4];
979 }
980 else if (skip == 3) {
981 l[2] = & _links[5];
982 }
983 else if (skip == 4) {
984 l[1] = & _links[2];
985 }
986 else if (skip == 5) {
987 l[2] = & _links[4];
988 }
989 else if (skip == 6) {
990 l[2] = & _links[3];
991 }
992 else if (skip == 7) {
993 l[1] = & _links[3];
994 l[2] = & _links[4];
995 }
996 else if (skip == 8) {
997 l[2] = & _links[5];
998 }
999 else if (skip == 9) {
1000 l[1] = & _links[4];
1001 }
1002 else if (skip == 10) {
1003 l[0] = & _links[1];
1004 }
1005 else if (skip == 11) {
1006 l[1] = & _links[3];
1007 }
1008 else if (skip == 12) {
1009 l[2] = & _links[4];
1010 }
1011 else if (skip == 13) {
1012 l[1] = & _links[2];
1013 }
1014 else if (skip == 14) {
1015 l[2] = & _links[3];
1016 }
1017 else if (skip == 15) {
1018 l[2] = & _links[5];
1019 }
1020 else if (skip == 16) {
1021 l[0] = & _links[2];
1022 l[1] = & _links[3];
1023 }
1024 else if (skip == 17) {
1025 l[2] = & _links[3];
1026 }
1027 else if (skip == 18) {
1028 l[1] = & _links[4];
1029 l[2] = & _links[5];
1030 }
1031 else if (skip == 19) {
1032 l[0] = & _links[3];
1033 }
1034
1035 unsigned n[3];
1036 for (unsigned i = 0; i < 3; i++)
1037 n[i] = l[i]->length();
1038
1039 for (unsigned i0 = 0; i0 < n[0]; i0++) {
1040 for (unsigned i1 = 0; i1 < n[1]; i1++) {
1041 for (unsigned i2 = 0; i2 < n[2]; i2++) {
1042 AList<TMLink> forLine;
1043 forLine.append((* l[0])[i0]);
1044 forLine.append((* l[1])[i1]);
1045 forLine.append((* l[2])[i2]);
1046
1047 TMLine & line = * new TMLine(forLine);
1048 int err = line.fit();
1049 if (err) {
1050 delete & line;
1051 continue;
1052 }
1053
1054 lines.append(line);
1055 }
1056 }
1057 }
1058 }
1059
1060#ifdef TRKRECO_DEBUG
1061 std::cout << " ... searchLines3 : # of lines found = "
1062 << lines.length() << std::endl;
1063#endif
1064
1065 return lines;
1066}
1067
1068AList<TMLine>
1069TBuilder::searchLines2(void) const {
1070 AList<TMLine> lines;
1071
1072 const AList<TMLink> * l[2];
1073 for (unsigned skip = 0; skip < 15; skip++) {
1074 if (skip == 0) {
1075 l[0] = & _links[0];
1076 l[1] = & _links[1];
1077 }
1078 else if (skip == 1) {
1079 l[1] = & _links[2];
1080 }
1081 else if (skip == 2) {
1082 l[1] = & _links[3];
1083 }
1084 else if (skip == 3) {
1085 l[1] = & _links[4];
1086 }
1087 else if (skip == 4) {
1088 l[1] = & _links[5];
1089 }
1090 else if (skip == 5) {
1091 l[0] = & _links[1];
1092 l[1] = & _links[2];
1093 }
1094 else if (skip == 6) {
1095 l[1] = & _links[3];
1096 }
1097 else if (skip == 7) {
1098 l[1] = & _links[4];
1099 }
1100 else if (skip == 8) {
1101 l[1] = & _links[5];
1102 }
1103 else if (skip == 9) {
1104 l[0] = & _links[2];
1105 l[1] = & _links[3];
1106 }
1107 else if (skip == 10) {
1108 l[1] = & _links[4];
1109 }
1110 else if (skip == 11) {
1111 l[1] = & _links[5];
1112 }
1113 else if (skip == 12) {
1114 l[0] = & _links[3];
1115 l[1] = & _links[4];
1116 }
1117 else if (skip == 13) {
1118 l[1] = & _links[5];
1119 }
1120 else if (skip == 14) {
1121 l[0] = & _links[4];
1122 }
1123
1124 unsigned n[2];
1125 for (unsigned i = 0; i < 2; i++)
1126 n[i] = l[i]->length();
1127
1128 for (unsigned i0 = 0; i0 < n[0]; i0++) {
1129 for (unsigned i1 = 0; i1 < n[1]; i1++) {
1130 AList<TMLink> forLine;
1131 forLine.append((* l[0])[i0]);
1132 forLine.append((* l[1])[i1]);
1133
1134 TMLine & line = * new TMLine(forLine);
1135 int err = line.fit();
1136 if (err) {
1137 delete & line;
1138 continue;
1139 }
1140
1141 lines.append(line);
1142 }
1143 }
1144 }
1145
1146#ifdef TRKRECO_DEBUG
1147 std::cout << " ... searchLines2 : # of lines found = "
1148 << lines.length() << std::endl;
1149#endif
1150
1151 return lines;
1152}
1153
1154AList<TMLine>
1155TBuilder::searchLines1(void) const {
1156 AList<TMLine> lines;
1157
1158 TMLine & line = * new TMLine(_forLine);
1159 int err = line.fit();
1160 if (err) {
1161 delete & line;
1162 }
1163 else {
1164 lines.append(line);
1165 }
1166
1167#ifdef TRKRECO_DEBUG
1168 std::cout << " ... searchLines1 : # of lines found = "
1169 << lines.length() << std::endl;
1170#endif
1171
1172 return lines;
1173}
1174
1175TMLine
1176TBuilder::searchLine(const TMLine & initialLine) const {
1177#ifdef TRKRECO_WINDOW
1178 std::string old = sz.text();
1179#endif
1180
1181 static float _szLinkMaxDistance = 5;
1182
1183 unsigned nLinks = _forLine.length();
1184 TMLine line0 = initialLine;
1185 line0.remove(initialLine.links());
1186#ifdef TRKRECO_DEBUG
1187 line0.fitted(true);
1188#endif
1189 unsigned nGoodLinks0 = 0;
1190 AList<TMLink> targets;
1191 unsigned nLoops = 0;
1192 while (1) {
1193 ++nLoops;
1194
1195 unsigned nGoodLinks = 0;
1196 for (unsigned i = 0; i < nLinks; i++) {
1197 float distance = line0.distance(* _forLine[i]);
1198 if (distance < _szLinkDistance)
1199 targets.append(_forLine[i]);
1200 if (distance < _szLinkMaxDistance) {
1201 ++nGoodLinks;
1202 }
1203 }
1204
1205#ifdef TRKRECO_DEBUG
1206 std::cout << " ... searchLine : # of close hits(last) = " << nGoodLinks0
1207 << ", # of close hits = " << nGoodLinks << std::endl;
1208#endif
1209
1210 //...Check condition...
1211 if (nGoodLinks0 == nGoodLinks) break;
1212
1213 //...Do fit...
1214 TMLine newLine(targets);
1215 static TRobustLineFitter fitter("robust fitter");
1216 int err = fitter.fit(newLine);
1217 if (err) {
1218#ifdef TRKRECO_WINDOW
1219 std::cout << " ... searchLine : failed to fit" << std::endl;
1220 std::string s = "line search failed";
1221 sz.text(s);
1222 sz.wait();
1223#endif
1224 break;
1225 }
1226
1227 //...To protect infinite loop...
1228 if (nLoops > 10) {
1229#ifdef TRKRECO_DEBUG_DETAIL
1230 std::cout << " reached to max # of loops(10) : break";
1231#endif
1232 break;
1233 }
1234
1235#ifdef TRKRECO_WINDOW
1236 std::string ss = ", # of close hits(last) = " + itostring(int(nGoodLinks0))
1237 + ", # of close hits = " + itostring(int(nGoodLinks));
1238 sz.append(line0, leda_brown);
1239 sz.append(newLine, leda_green);
1240 sz.text(old + ", nLoops = " + itostring(int(nLoops)) + ss);
1241 sz.wait();
1242 sz.remove(line0);
1243 sz.remove(newLine);
1244#endif
1245
1246 line0 = newLine;
1247 nGoodLinks0 = nGoodLinks;
1248 targets.removeAll();
1249 }
1250
1251 return line0;
1252}
1253
1254TTrack *
1255TBuilder::build(TTrack & t, const TMLine & l) const {
1256 AList<TMLink> links = l.links();
1257 AList<TMLink> toRemove;
1258 unsigned n = links.length();
1259 for (unsigned i = 1; i < n; i++) {
1260 if (links[i - 1]->link() == links[i]->link()) {
1261 toRemove.append(links[i]);
1262 continue;
1263 }
1264 if (i < 2)
1265 continue;
1266 if (links[i - 2]->link() == links[i]->link())
1267 toRemove.append(links[i]);
1268 }
1269 links.remove(toRemove);
1270
1271 //...Pick up links...
1272 n = links.length();
1273 for (unsigned i = 0; i < n; i++)
1274 t.append(* links[i]->link());
1275
1276 Vector a(5);
1277 a = t.helix().a();
1278 a[3] = l.b();
1279 //a[4] = t.charge() * l.a();//yzhang delete 2012-05-08
1280 // //yzhang add 2012-05-08
1281 if(_fitter.getMagneticFieldPointer() == NULL){
1282 std::cout<<" "<<__FILE__<<" "<<__LINE__<<" Could not get MagneticFieldSvc "<<std::endl;
1283 }
1284 const double Bz = -1000 * _fitter.getMagneticFieldPointer()->getReferField();
1285 if(Bz>0){
1286 a[4] = t.charge() * l.a();//yzhang add 2012-05-08 ???
1287 }else{
1288 a[4] = -1. * t.charge() * l.a();//yzhang add 2012-05-08 ???
1289 }
1290 //zhangy
1291 t._helix->a(a);
1292
1293//#ifdef TRKRECO_WINDOW
1294// TTrack * td = new TTrack(t);
1295// sz.append(* td, leda_green);
1296//#endif
1297#ifdef TRKRECO_DEBUG_DETAIL
1298 static unsigned nTrk = 0;
1299 ++nTrk;
1300// for (unsigned i = 0; i < t.nCores(); i++)
1301// if (t.cores()[i]->wire()->name() == "43=220") {
1302// std::cout << "43=220 removed" << std::endl;
1303// t.remove(* t.cores()[i]);
1304// }
1305 HepVector3D pp;
1306 if (links[0]->hit()->mc())
1307 pp = links[0]->hit()->mc()->hep()->p().vect();
1308 const HepVector3D p0 = pp;
1309// if (nTrk == 1) {
1310// Helix tmph(Point3D(0, 0, 0), p0, +1);
1311// * t._helix = tmph;
1312// }
1313 t.dump("detail", "before fit");
1314 std::cout << "Pdif mag=" << (t.p() - p0).mag() << std::endl;
1315// t.links().removeAll();
1316// for (unsigned i = 0; i < BsCouTab(RECMDC_WIRHIT
1317// }
1318#endif
1319
1320 //...Refine...
1321 AList<TMLink> bad;
1322 int err = _fitter.fit(t);
1323
1324#ifdef TRKRECO_DEBUG_DETAIL
1325 t.dump("detail", "after fit");
1326 std::cout << "Pdif mag=" << (t.p() - p0).mag() << std::endl;
1327#endif
1328
1329 t.refine(bad, _maxSigmaStereo * 100.);
1330 err = _fitter.fit(t);
1331 t.refine(bad, _maxSigmaStereo * 10.);
1332 err = _fitter.fit(t);
1333 t.refine(bad, _maxSigmaStereo);
1334#ifdef TRKRECO_DEBUG_DETAIL
1335// if (nTrk == 2) {
1336// Helix tmph(Point3D(0, 0, 0), p0, +1);
1337// * t._helix = tmph;
1338// std::cout << "initial mom" << p0 << " is set" << std::endl;
1339// }
1340#endif
1341 err = _fitter.fit(t);
1342
1343#ifdef TRKRECO_DEBUG_DETAIL
1344 t.dump("detail"," ");
1345 std::cout << "Pdif mag=" << (t.p() - p0).mag() << std::endl;
1346#endif
1347
1348 //...Termination...
1349 return & t;
1350}
1351
1352bool
1353TBuilder::initializeForStereo(const TTrack & t,
1354 const AList<TSegment> & segments,
1355 const AList<TSegment> & badSegments) const {
1356 _nSuperLayers = 0;
1357 for (unsigned i = 0; i < 5; i++)
1358 _nHits[i] = 0;
1359 //...sz position of segments...
1360 unsigned nSegments = segments.length();
1361// cout<<" used Segments in stereo = "<<nSegments<<endl;
1362 for (unsigned i = 0; i < nSegments; i++) {
1363 TMLink & l = * new TMLink();
1364 int err = t.szPosition(* segments[i], l);
1365 //...Remove if sz can not be calculcated...
1366 if (err) {
1367 delete & l;
1368 continue;
1369 }
1370// _links[l.wire()->superLayerId() / 2].append(l);
1371 if (!_links[l.wire()->axialStereoLayerId()/4].hasMember(l))
1372 _links[l.wire()->axialStereoLayerId()/4].append(l);
1373//? _allLinks.append(segments[i]->links());
1374// _allLinks.append(segments[i]->cores());
1375 for (unsigned j = 0; j < segments[i]->cores().length(); ++j) {
1376 TMLink & cL = * segments[i]->cores()[j];
1377 if(!_allLinks.hasMember(cL)) _allLinks.append(cL);
1378 }
1379 }
1380
1381 //...Count a number of super layers...
1382 for (unsigned i = 0; i < 6; i++) {
1383 if (_links[i].length() > 0) {
1384 ++_nSuperLayers;
1385 }
1386// else {
1387// _allLinks.append(badSegments[i]->links());
1388// #ifdef TRKRECO_DEBUG_DETAIL
1389// std::cout << " ... stereo super layer " << i * 2 + 1
1390// << " has no link,"
1391// << badSegments[i]->links().length() << " (bad) links added"
1392// << std::endl;
1393// badSegments[i]->dump("hits sort flag", " ");
1394// #endif
1395// }
1396 }
1397
1398 //...Append bad links also...
1399 if (badSegments.length()) {
1400 for (unsigned i = 0; i < 6; i++) {
1401 if (badSegments[i]->links().length()) {
1402 _allLinks.append(badSegments[i]->links());
1403#ifdef TRKRECO_DEBUG_DETAIL
1404 std::cout << " ... bad links added for stereo super layer "
1405 << i * 2 + 1 << std::endl;
1406 badSegments[i]->dump("hits sort flag", " ");
1407#endif
1408 }
1409 }
1410 }
1411
1412 //...sz position of links...
1413 unsigned nCores = _allLinks.length();
1414 if (nCores < _minNCores) {
1415#ifdef TRKRECO_DEBUG_DETAIL
1416 std::cout << " ... initializeForStereo : # of cores(=" << nCores
1417 << ") is less then " << _minNCores << std::endl;
1418#endif
1419 return false;
1420 }
1421 for (unsigned j = 0; j < nCores; j++) {
1422 TMLink & link = * _allLinks[j];
1423 for (unsigned i = 0; i < 2; i++) {
1424 TMLink & tt = * new TMLink(link);
1425 unsigned stateR=(link.hit()->state())&WireHitPatternRight;
1426 unsigned stateL=(link.hit()->state())&WireHitPatternLeft;
1427 bool boolR=(stateR==WireHitPatternRight);
1428 bool boolL=(stateL==WireHitPatternLeft);
1429// std::cout<<"zsl TBuilder::initializeForStereo::boolR:"<<boolR<<", boolL:"<<boolL<<endl;
1430
1431 tt.leftRight(i);
1432 int err = t.szPosition(tt);
1433// std::cout<<"zsl TBuilder::initializeForStereo::err:"<<err<<endl;
1434// std::cout<<"----------------------"<<endl;
1435
1436 if (err) {
1437 delete & tt;
1438 continue;
1439 }
1440 tt.link(& link);
1441 _forLine.append(tt);
1442 }
1443 }
1444
1445 //...Count # of axial hits...
1446 unsigned nA = t.cores().length();
1447 for (unsigned i = 0; i < nA; i++)
1448// ++_nHits[t.cores()[i]->wire()->superLayerId() / 2];
1449 ++_nHits[t.cores()[i]->wire()->axialStereoLayerId()/4];
1450
1451#ifdef TRKRECO_DEBUG
1452 std::cout << " ... initializeForStereo : axial super layer usage = ";
1453 for (unsigned i = 0; i < 5; i++)
1454 std::cout << _nHits[i] << ",";
1455 std::cout << std::endl
1456 << " : # of stereo super layers="
1457 << _nSuperLayers << std::endl;
1458#endif
1459
1460 return true;
1461}
1462
1463unsigned
1464TBuilder::stereoQuality(const AList<TMLink> & links) const {
1465 unsigned n[6] = {0, 0, 0, 0, 0, 0};
1466
1467 //...Check superlayers...
1468 unsigned nLinks = links.length();
1469 for (unsigned i = 0; i < nLinks; i++) {
1470 const TMLink & l = * links[i];
1471 if (l.wire()->axial()) continue;
1472// ++n[l.wire()->superLayerId() / 2];
1473// cout<<l.wire()->axialStereoLayerId()<<endl;
1474 ++n[l.wire()->axialStereoLayerId()/4];
1475 }
1476
1477 unsigned nTotal = 0;
1478 unsigned nMissing = 0;
1479 unsigned innermostStereo = 5;
1480 unsigned outermostStereo = 0;
1481 unsigned innermostAxial = 4;
1482 unsigned outermostAxial = 0;
1483 unsigned innerSL = 6;
1484 unsigned outerSL = 0;
1485 for (unsigned i = 0; i < 6; i++) {
1486 if (_links[i].length() > 1) {
1487 if(i > outermostStereo) outermostStereo = i;
1488 if(i < innermostStereo) innermostStereo = i;
1489 }
1490 if (n[i] > 1) ++nTotal;
1491 }
1492 for (unsigned i = 0; i < 5; i++) {
1493 if (_nHits[i] > 1) {
1494 if (i > outermostAxial) outermostAxial = i;
1495 if (i < innermostAxial) innermostAxial = i;
1496 }
1497 }
1498
1499 if(innermostStereo > 1 && innermostAxial < 3) innerSL = 2;
1500 else innerSL = innermostStereo;
1501 if(outermostStereo > 1 && outermostAxial > 2) outerSL = 5;
1502 else outerSL = outermostStereo;
1503
1504 for (unsigned i = 0; i < 6; i++) {
1505 if (_links[i].length() > 0) {
1506//zsl 051014 if ((_nHits[i] > 1) && (_nHits[i + 1] > 1))
1507 if (i <= outerSL && i >= innerSL)
1508 if (n[i] < 2)
1509 ++nMissing;
1510 }
1511 if (n[i] > 1) ++nTotal;
1512 }
1513 unsigned toBeReturn = 0;
1514 if (nMissing <= 1) toBeReturn = 2;
1515 else if (nMissing == 2) toBeReturn = 1;
1516 if (nTotal < 2) toBeReturn = 0;
1517
1518#ifdef TRKRECO_DEBUG
1519 std::cout << " ... stereoQuality : axial ";
1520 for (unsigned i = 0; i < 6; i++)
1521 std::cout << _nHits[i] << ",";
1522 std::cout << std::endl
1523 << " : stereo ";
1524 for (unsigned i = 0; i < 5; i++)
1525 std::cout << n[i] << ",";
1526 std::cout << " : total " << nTotal << " super layers, "
1527 << " : quality=" << toBeReturn << std::endl;
1528#endif
1529
1530 return toBeReturn;
1531}
1532
1533
1534TTrack *
1535TBuilder::buildStereoNew(const TTrack & t,
1536 AList<TSegment> & segments,
1537 AList<TSegment> & badSegments) const {
1538
1539#ifdef TRKRECO_DEBUG
1540 std::cout << "... building stereo by segments(new) : # of segments = ";
1541 std::cout << segments.length() << std::endl;
1542 for (unsigned i = 0; i < segments.length(); i++)
1543 segments[i]->dump("hits sort flag", " ");
1544#endif
1545#ifdef TRKRECO_WINDOW
1546 sz.clear();
1547#endif
1548 TTrack * bestCandidate = NULL;
1549 AList< AList<TMLink> > poorSeeds;
1550 bool ok = initializeForStereo(t, segments, badSegments);
1551 unsigned nSuperLayers = _nSuperLayers + 1;
1552// std::cout<<"nSuperLayers in buildStereoNew:"<<_nSuperLayers<<endl;
1553// cout<<"ok of initializeForStereo: "<<ok<<endl;
1554
1555 if (! ok) goto endOfBuilding;
1556 //...Main loop...
1557 while (--nSuperLayers) {
1558
1559 //...Initial line search...
1560 AList<TMLine> initialLines = searchInitialLines(nSuperLayers);
1561#ifdef TRKRECO_WINDOW
1562 sz.clear();
1563 sz.skip(false);
1564 sz.mode(2);
1565 sz.appendSz(t, segments, leda_black);
1566 AList<TSegment> tmps = segments;
1567 std::string s = "# of segments = " + itostring(int(segments.length()));
1568 sz.appendSz(t, _allLinks, leda_black);
1569 sz.append(_forLine, leda_brown);
1570 s = "nSprLyr=" + itostring(int(nSuperLayers)) +
1571 ", # of initial lines = " + itostring(int(initialLines.length()));
1572 for (unsigned i = 0; i < initialLines.length(); i++)
1573 sz.append(* initialLines[i], leda_red);
1574 sz.text(s);
1575 sz.wait();
1576#endif
1577// cout<<"length of initialLines = "<<initialLines.length()<<endl;
1578 if (initialLines.length() == 0) continue;
1579 //...Line loop...
1580 bool found = false;
1581 unsigned nInitialLines = initialLines.length();
1582// cout<<"nInitialLines:"<<nInitialLines<<endl;
1583
1584 for (unsigned i = 0; i < nInitialLines; i++) {
1585
1586 //...Linear fit...
1587 const TMLine & line = * initialLines[i];
1588 TMLine newLine = searchLine(line);
1589
1590 //...Skip if this is a seed of the poor result...
1591 if (poorSeeds.length()) {
1592 bool poorCase = false;
1593 for (unsigned j = 0; j < poorSeeds.length(); j++) {
1594 if (poorSeeds[j]->length() == newLine.links().length()) {
1595 AList<TMLink> tmp = * poorSeeds[j];
1596 tmp.remove(newLine.links());
1597 if (tmp.length() == 0) {
1598#ifdef TRKRECO_DEBUG_DETAIL
1599 std::cout << " ... This is a poor seed :"
1600 << " skipped"
1601 << " : # of poor seeds = "
1602 << poorSeeds.length() << std::endl;
1603#endif
1604 poorCase = true;
1605 break;
1606 }
1607 }
1608 }
1609 if (poorCase) continue;
1610 }
1611
1612 //...Is this a good line...
1613 if (! stereoQuality(newLine.links()))
1614 continue;
1615
1616 //...3D fit...
1617 TTrack * t3d = new TTrack(t);
1618 t3d = build(* t3d, newLine);
1619 if (t3d == NULL) continue;
1620
1621 //...Check quality...
1622 unsigned quality = stereoQuality(t3d->links());
1623// cout<<"point3, quality:"<<quality<<endl;
1624
1625 //...Best case...
1626 if (quality == 2) {
1627#ifdef TRKRECO_WINDOW
1628 sz.text("stereo finished");
1629 sz.oneShot(* t3d, leda_blue);
1630#endif
1631 if (bestCandidate) delete bestCandidate;
1632 bestCandidate = t3d;
1633 found = true;
1634 break;
1635 }
1636
1637 //...Poor case...
1638 AList<TMLink> * tmpL = new AList<TMLink>();
1639 tmpL->append(newLine.links());
1640 poorSeeds.append(tmpL);
1641 if (quality == 0) {
1642#ifdef TRKRECO_WINDOW
1643 sz.text("this candidate discarded");
1644 sz.oneShot(* t3d, leda_black);
1645#endif
1646 delete t3d;
1647 continue;
1648 }
1649
1650 //...Not enough...
1651 if (bestCandidate) {
1652 if (bestCandidate->cores().length() <
1653 t3d->cores().length()) {
1654#ifdef TRKRECO_WINDOW
1655 sz.text("new candidate");
1656 sz.append(* bestCandidate, leda_brown);
1657 sz.oneShot(* t3d, leda_green);
1658 sz.remove(* bestCandidate);
1659#endif
1660 delete bestCandidate;
1661 bestCandidate = t3d;
1662 }
1663 else {
1664#ifdef TRKRECO_WINDOW
1665 sz.text("this candidate discarded");
1666 sz.oneShot(* t3d, leda_black);
1667#endif
1668 delete t3d;
1669 }
1670 }
1671 else {
1672 bestCandidate = t3d;
1673#ifdef TRKRECO_WINDOW
1674 sz.text("new candidate");
1675 sz.oneShot(* bestCandidate, leda_green);
1676 sz.remove(* t3d);
1677#endif
1678 }
1679 }
1680
1681 //...Termination of a loop...
1682 HepAListDeleteAll(initialLines);
1683 if (found) break;
1684 }
1685
1686endOfBuilding:
1687 _allLinks.removeAll();
1688 for (unsigned i = 0; i < 6; i++)
1689 HepAListDeleteAll(_links[i]);
1690 HepAListDeleteAll(_forLine);
1691
1692#ifdef TRKRECO_DEBUG_DETAIL
1693 std::cout << " ... # of poor seeds = " << poorSeeds.length() << std::endl;
1694#endif
1695#ifdef TRKRECO_WINDOW
1696 if (bestCandidate == NULL) {
1697 sz.text("3D failed");
1698 sz.wait();
1699 }
1700#endif
1701
1702 //...Check used segments...
1703 if (bestCandidate) {
1704 AList<TMLink> usedLinks = StereoHits(bestCandidate->links());
1705 for (unsigned i = 0; i < segments.length(); i++) {
1706 TSegment & segment = * segments[i];
1707 AList<TMLink> used = Links(segment, * bestCandidate);
1708 if (used.length()) {
1709 bestCandidate->segments().append(segment);
1710 segment.tracks().append(bestCandidate);
1711 }
1712 }
1713 }
1714
1715 if (poorSeeds.length())
1716 HepAListDeleteAll(poorSeeds);
1717
1718#ifdef TRASA_DEBUG_DETAIL
1719 if (bestCandidate == NULL)
1720 std::cout << "... building stereo(new) failed" << std::endl;
1721 else
1722 std::cout << "... building stereo(new) ok" << std::endl;
1723#endif
1724 return bestCandidate;
1725}
n
const Int_t n
Definition DataBase/tau_mode.c:65
s
XmlRpcServer s
Definition HelloServer.cpp:11
WireHitPatternLeft
#define WireHitPatternLeft
Definition TMDCWireHit.h:33
WireHitPatternRight
#define WireHitPatternRight
Definition TMDCWireHit.h:34
SortByB
int SortByB(const void *av, const void *bv)
Sorter.
Definition TMLine.cxx:648
Links
AList< TMLink > Links(const TSegment &s, const TTrack &t)
returns AList of TMLink used for a track.
Definition TSegment.cxx:963
NULL
#define NULL
Definition VertexExtrapolate.h:14
t
TTree * t
Definition binning.cxx:23
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:263
IMagneticFieldSvc::getReferField
virtual double getReferField()=0
TBuilder::initialLine
TMLine * initialLine(const TTrack &, AList< TSegment > &) const
makes a line.
Definition TBuilder.cxx:315
TBuilder::searchLines5
AList< TMLine > searchLines5(void) const
Definition TBuilder.cxx:798
TBuilder::searchLines1
AList< TMLine > searchLines1(void) const
Definition TBuilder.cxx:1155
TBuilder::initialLineOld
TMLine * initialLineOld(const TTrack &, AList< TSegment > &) const
Definition TBuilder.cxx:405
TBuilder::removeFarSegment
void removeFarSegment(const TMLine &, AList< TSegment > &, AList< TMLink > &) const
Definition TBuilder.cxx:510
TBuilder::searchLines6
AList< TMLine > searchLines6(void) const
Definition TBuilder.cxx:753
TBuilder::buildRphi
TTrack * buildRphi(const AList< TMLink > &) const
builds a r/phi track from TMLinks or from Segments.
Definition TBuilder.cxx:107
TBuilder::initialLine2
TMLine * initialLine2(const TTrack &, const AList< TMLink > &) const
Definition TBuilder.cxx:296
TBuilder::~TBuilder
virtual ~TBuilder()
Destructor.
Definition TBuilder.cxx:58
TBuilder::build
TTrack * build(TTrack &t, const TMLine &l) const
Definition TBuilder.cxx:1255
TBuilder::salvage
void salvage(TTrack &t, AList< TMLink > &hits) const
salvages hits.
Definition TBuilder.cxx:211
TBuilder::searchLines2
AList< TMLine > searchLines2(void) const
Definition TBuilder.cxx:1069
TBuilder::TBuilder
TBuilder(const std::string &name, float maxSigma, float maxSigmaStereo, float salvageLevel, float szSegmentDistance, float szLinkDistance, unsigned fittingFlag)
Constructor with salvage level.
Definition TBuilder.cxx:36
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
TBuilder::selectStereoSegment
AList< TSegment > selectStereoSegment(const TMLine &line, const AList< TSegment > &list, const AList< TMLink > &szList) const
Definition TBuilder.cxx:497
TBuilder::initialLine1
TMLine * initialLine1(const TTrack &, const AList< TSegment > &, const AList< TMLink > &) const
Definition TBuilder.cxx:247
TBuilder::dump
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TBuilder.cxx:62
TBuilder::searchLines4
AList< TMLine > searchLines4(void) const
Definition TBuilder.cxx:867
TBuilder::searchInitialLines
AList< TMLine > searchInitialLines(unsigned nSuperLayers) const
Definition TBuilder.cxx:724
TBuilder::searchLines3
AList< TMLine > searchLines3(void) const
Definition TBuilder.cxx:964
TBuilder::searchLine
TMLine searchLine(const TMLine &initialLine) const
Definition TBuilder.cxx:1176
TCircle
A class to represent a circle in tracking.
Definition TCircle.h:42
TCircle::dump
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TCircle.cxx:33
THelixFitter::freeT0
bool freeT0(void) const
sets/returns free T0 flag.
Definition THelixFitter.h:275
THelixFitter::tof
bool tof(void) const
sets/returns propagation-delay correction flag.
Definition THelixFitter.h:207
THelixFitter::fit
int fit(TTrackBase &) const
Definition THelixFitter.h:231
THelixFitter::sag
bool sag(void) const
sets/returns sag correction flag.
Definition THelixFitter.h:181
THelixFitter::getMagneticFieldPointer
IMagneticFieldSvc * getMagneticFieldPointer(void) const
Definition THelixFitter.h:98
THelixFitter::propagation
bool propagation(void) const
sets/returns propagation-delay correction flag.
Definition THelixFitter.h:193
TMDCWire::axial
bool axial(void) const
returns true if this wire is in an axial layer.
Definition TMDCWire.h:348
TMDCWire::axialStereoLayerId
unsigned axialStereoLayerId(void) const
returns id of axial or stereo id.
Definition TMDCWire.h:360
TMDCWire::forwardPosition
const HepPoint3D & forwardPosition(void) const
returns position in forward endplate.
Definition TMDCWire.h:300
TMDCWire::backwardPosition
const HepPoint3D & backwardPosition(void) const
returns position in backward endplate.
Definition TMDCWire.h:306
TMLine
A class to represent a track in tracking.
Definition TMLine.h:40
TMLine::distance
double distance(const TMLink &) const
returns distance to a position of TMLink itself. (not to a wire)
Definition TMLine.h:165
TMLine::a
double a(void) const
returns coefficient a.
Definition TMLine.h:147
TMLine::b
double b(void) const
returns coefficient b.
Definition TMLine.h:156
TRobustLineFitter
A class to fit a TTrackBase object to a line.
Definition TRobustLineFitter.h:28
TRobustLineFitter::fit
virtual int fit(TTrackBase &) const
Definition TRobustLineFitter.cxx:45
TSegment
A class to relate TMDCWireHit and TTrack objects.
Definition TSegment.h:43
TSegment::tracks
AList< TTrack > & tracks(void)
Definition TSegment.h:360
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::refine
virtual void refine(AList< TMLink > &list, double maxSigma)
removes bad points by pull. The bad points are removed from the track, and are returned in 'list'.
Definition TTrackBase.cxx:170
TTrackBase::fitted
bool fitted(void) const
returns true if fitted.
Definition TTrackBase.h:222
TTrackBase::append
void append(TMLink &)
appends a TMLink.
Definition TTrackBase.cxx:362
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::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
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::helix
const Helix & helix(void) const
returns helix parameter.
Definition TTrack.h:477
TTrack::szPosition
int szPosition(TMLink &link) const
calculates arc length and z for a stereo hit.
Definition TTrack.cxx:3372
TTrack::charge
double charge(void) const
returns charge.
Definition TTrack.h:504