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