BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
TBuilder0.cxx
Go to the documentation of this file.
1//-----------------------------------------------------------------------------
2// $Id: TBuilder0.cxx,v 1.9 2010/03/31 09:58:59 liucy Exp $
3//-----------------------------------------------------------------------------
4// Filename : TBuilder0.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#ifdef TRKRECO_DEBUG_DETAIL
14# ifndef TRKRECO_DEBUG
15# define TRKRECO_DEBUG
16# endif
17#endif
18#include "TrkReco/TBuilder0.h"
19#include "TrkReco/TCircle.h"
20#include "TrkReco/TLine0.h"
21#include "TrkReco/TMDC.h"
22#include "TrkReco/TMDCWire.h"
23#include "TrkReco/TMDCWireHit.h"
24#include "TrkReco/TMLink.h"
25#include "TrkReco/TTrack.h"
26////#include "TrkReco/TSegment0.h"
27#ifdef TRKRECO_WINDOW
28# include "TrkReco/TWindow.h"
29#endif
30
31TBuilder0::TBuilder0( const std::string& a )
32 : _name( a )
33 , _fitter( "TBuilder0 Fitter" )
34 , _salvageLevel( 30. )
35 , _stereoZ3( 20. )
36 , _stereoZ4( 20. )
37 , _stereoChisq3( 15. )
38 , _stereoChisq4( 9. )
39 , _stereoMaxSigma( 30. ) {}
40
41TBuilder0::TBuilder0( const std::string& a, float salvageLevel )
42 : _name( a )
43 , _fitter( "TBuilder0 Fitter" )
44 , _salvageLevel( salvageLevel )
45 , _stereoZ3( 20. )
46 , _stereoZ4( 20. )
47 , _stereoChisq3( 15. )
48 , _stereoChisq4( 9. )
49 , _stereoMaxSigma( 30. ) {}
50
51TBuilder0::TBuilder0( const std::string& a, float stereoZ3, float stereoZ4, float stereoChisq3,
52 float stereoChisq4, float stereoMaxSigma, unsigned corrections,
53 float salvageLevel )
54 : _name( a )
55 , _fitter( "TBuilder0 Fitter" )
56 , _stereoZ3( stereoZ3 )
57 , _stereoZ4( stereoZ4 )
58 , _stereoChisq3( stereoChisq3 )
59 , _stereoChisq4( stereoChisq4 )
60 , _stereoMaxSigma( stereoMaxSigma )
61 , _salvageLevel( salvageLevel ) {
62 // if (corrections & 1) _fitter.sag(true);
63 // if (corrections & 2) _fitter.propagation(true);
64 // if (corrections & 4) _fitter.tof(true);
65 // if (corrections & 8) _fitter.freeT0(true);
66}
67
68TBuilder0::~TBuilder0() {}
69
70void TBuilder0::dump( const std::string& msg, const std::string& pre ) const {}
71
72TTrack* TBuilder0::buildRphi( const AList<TMLink>& list ) const {
73#ifdef TRKRECO_DEBUG_DETAIL
74 std::cout << _name << " ... building a rphi track" << std::endl;
75 std::cout << _name << " ... selecting good hits" << std::endl;
76// TTrackBase tmp;
77// tmp.append(list);
78// rphiWindow->append(tmp, leda_red);
79#endif
80
81 //...Check # of links...
82 if ( list.length() < _circleSelector.nLinks() )
83 {
84#ifdef TRKRECO_DEBUG_DETAIL
85 std::cout << _name << " ... rejected by nLinks(";
86 std::cout << list.length() << ") < ";
87 std::cout << _circleSelector.nLinks() << std::endl;
88#endif
89 }
90
91 //...Select core hits...
92 AList<TMLink> cores = list;
93 selectHits( cores );
94 if ( cores.length() < 5 ) cores = list;
95 // cout<<"list: "<<list.length()<<" core: "<<cores.length()<<endl;
96
97 //...Core check...
98#ifdef TRKRECO_DEBUG_DETAIL
99 std::cout << _name << " ... checking cores : cores=" << std::endl;
100 Dump( cores, "detail" );
101#endif
102 unsigned sLinks = SuperLayer( list );
103 unsigned sUsed = SuperLayer( cores );
104#ifdef TRKRECO_DEBUG_DETAIL
105 std::cout << " super layer ptn=" << sLinks;
106 std::cout << ",super layer used=" << sUsed << std::endl;
107#endif
108 if ( sLinks != sUsed )
109 {
110#ifdef TRKRECO_DEBUG_DETAIL
111 std::cout << _name << " ... appending hits to cores" << std::endl;
112#endif
113 unsigned diff = sLinks - sUsed;
114 unsigned asl[5] = { 2, 3, 4, 9, 10 };
115 for ( unsigned j = 0; j < 5; j++ )
116 {
117 // if (diff & (1 << ((5 - j) * 2))) {
118 // if ( diff & ( 1 << asl[5 - j] ) )
119 if ( diff & ( 1 << asl[4 - j] ) ) // yzhang debug for 720 FIXME
120 {
121#ifdef TRKRECO_DEBUG_DETAIL
122 std::cout << " super layer " << ( 5 - j ) * 2 << "searching";
123 std::cout << std::endl;
124#endif
125 // AList<TMLink> links = SameSuperLayer( list, asl[5 - j] );
126 AList<TMLink> links = SameSuperLayer( list, asl[4 - j] ); // yzhang debug for 720 FIXME
127 TMLink* best = NULL;
128 double bestD = 999.;
129 for ( unsigned i = 0; i < links.length(); i++ )
130 {
131 double dist = links[i]->hit()->drift();
132 if ( dist < 0.02 )
133 {
134#ifdef TRKRECO_DEBUG_DETAIL
135 std::cout << " " << links[i]->wire()->name();
136 std::cout << " appended (small dist)" << std::endl;
137#endif
138 cores.append( links[i] );
139 continue;
140 }
141 if ( dist < bestD )
142 {
143 best = links[i];
144 bestD = dist;
145 }
146 }
147 if ( best )
148 {
149 cores.append( best );
150#ifdef TRKRECO_DEBUG_DETAIL
151 std::cout << " " << best->wire()->name();
152 std::cout << " appended (best)" << std::endl;
153#endif
154 }
155 }
156 }
157 }
158
159 //...Check cores again...
160 unsigned nCores = cores.length();
161 AList<TMLink> realCores;
162 for ( unsigned i = 0; i < nCores; i++ )
163 {
164 TMLink* l = cores[i];
165 if ( l->hit()->state() & WireHitFittingValid ) realCores.append( l );
166 }
167 if ( NSuperLayers( realCores ) < 2 )
168 {
169#ifdef TRKRECO_DEBUG_DETAIL
170 std::cout << " ... rejected by small number of super layers" << std::endl;
171#endif
172 return NULL;
173 }
174 if ( NLayers( realCores ) < 5 )
175 {
176#ifdef TRKRECO_DEBUG_DETAIL
177 std::cout << " ... rejected by small number of layers" << std::endl;
178#endif
179 return NULL;
180 }
181
182 //...Make a circle...
183#ifdef TRKRECO_DEBUG_DETAIL
184 std::cout << _name << " ... making a circle : #cores=" << cores.length();
185 std::cout << std::endl;
186#endif
187 AList<TMLink> hits = list;
188 hits.remove( cores );
189 TCircle c( cores );
190
191 //...Test it...
192 if ( !_circleSelector.preSelect( c ) ) return NULL;
193
194 //...Fitting...
195 int err = c.fit();
196 if ( err < 0 )
197 {
198#ifdef TRKRECO_DEBUG_DETAIL
199 std::cout << " ... rejected by failure of the 1st fit : ";
200 std::cout << "err = " << err << std::endl;
201#endif
202 return NULL;
203 }
204
205 //...Test it...
206 if ( !_circleSelector.select( c ) ) return NULL;
207
208 //...Make a track...
209#ifdef TRKRECO_DEBUG_DETAIL
210 std::cout << _name << " ... making a track" << std::endl;
211#endif
212 TTrack* t = new TTrack( c );
213 if ( !_trackSelector.preSelect( *t ) )
214 {
215 delete t;
216 return NULL;
217 }
218
219 //...Fitting...
220 AList<TMLink> bad;
221 // cout<<"track0: "<<t->nLinks()<<endl;
222 err = fit( *t );
223 if ( err < 0 ) goto discard;
224 t->refine( bad, _trackSelector.maxSigma() * 100. );
225 // cout<<"track1: "<<t->nLinks()<<endl;
226 err = fit( *t );
227 if ( err < 0 ) goto discard;
228#ifdef TRKRECO_DEBUG_DETAIL
229 t->dump( "detail", " 1st> " );
230#endif
231
232 //...Test it...
233 if ( !_trackSelector.select( *t ) ) goto discard;
234
235 //...Try to append non-core hits...
236#ifdef TRKRECO_DEBUG_DETAIL
237 std::cout << _name << " ... appending non-core hits" << std::endl;
238#endif
239 t->appendByApproach( hits, sqrt( _trackSelector.maxSigma() ) );
240 err = fit( *t );
241 if ( err < 0 ) goto discard;
242 // cout<<" 2D track1: "<<t->links().length()<<endl;
243 t->refine( bad, _trackSelector.maxSigma() * 10. );
244 err = fit( *t );
245 if ( err < 0 ) goto discard;
246 // cout<<" 2D track2: "<<t->links().length()<<endl;
247 t->refine( bad, _trackSelector.maxSigma() );
248 err = fit( *t );
249 if ( err < 0 ) goto discard;
250// cout<<" 2D track3: "<<t->links().length()<<endl;
251#ifdef TRKRECO_DEBUG_DETAIL
252 t->dump( "detail", " 2nd> " );
253#endif
254
255 //...Test it...
256 if ( !_trackSelector.select( *t ) ) goto discard;
257
258 //...OK...
259#ifdef TRKRECO_DEBUG_DETAIL
260// rphiWindow->append(* t, leda_blue);
261// rphiWindow->draw();
262// rphiWindow->wait();
263// rphiWindow->remove(tmp);
264// rphiWindow->remove(* t);
265#endif
266
267 return t;
268
269 //...Something happened...
270discard:
271#ifdef TRKRECO_DEBUG_DETAIL
272 std::cout << " ... rejected by fitting failure : ";
273 std::cout << " err = " << err << std::endl;
274// rphiWindow->append(* t, leda_blue);
275// rphiWindow->draw();
276// rphiWindow->wait();
277// rphiWindow->remove(tmp);
278// rphiWindow->remove(* t);
279#endif
280 delete t;
281 return NULL;
282}
283
284void TBuilder0::selectHits( AList<TMLink>& list ) const {
285 AList<TMLink> bad;
286 for ( unsigned i = 0; i < list.length(); i++ )
287 {
288 unsigned state = list[i]->hit()->state();
289 if ( ( !( state & WireHitContinuous ) ) || ( !( state & WireHitIsolated ) ) )
290 {
291 bad.append( list[i] );
292 continue;
293 }
294/* if ((! (state & WireHitPatternLeft)) &&
295 (! (state & WireHitPatternRight))) {
296 bad.append(list[i]);
297 continue;
298 }
299*/ }
300
301#ifdef TRKRECO_DEBUG_DETAIL
302std::cout << _name << "(TBuilder0::selectHits) ... dump of candidates" << std::endl;
303Dump( list, "detail" );
304#endif
305
306list.remove( bad );
307}
308
309TTrack* TBuilder0::buildStereo0( TTrack& track, const AList<TMLink>& list ) const {
310#ifdef TRKRECO_DEBUG_DETAIL
311 std::cout << _name << "(stereo) ... dump of stereo candidate hits" << std::endl;
312 Dump( list, "sort flag", " " );
313#endif
314
315 //...Check # of links...
316 if ( list.length() < _lineSelector.nLinksStereo() )
317 {
318#ifdef TRKRECO_DEBUG_DETAIL
319 std::cout << _name << "(stereo) ... rejected by nLinks(";
320 std::cout << list.length() << ") < ";
321 std::cout << _lineSelector.nLinks() << std::endl;
322#endif
323 return NULL;
324 }
325
326 //...Calculate s and z for every links...
327 unsigned n = list.length();
328 AList<TMLink> forLine;
329 for ( unsigned i = 0; i < n; i++ )
330 {
331 TMLink* l = list[i];
332 TMLink* t = new TMLink( *l );
333
334 //...Assuming wire position...
335 t->leftRight( 2 );
336 int err = track.szPosition( *t );
337 if ( err )
338 {
339 delete t;
340 continue;
341 }
342
343 //...Store the sz link...
344 t->link( l );
345 forLine.append( t );
346 }
347
348#ifdef TRKRECO_DEBUG_DETAIL
349 std::cout << _name << "(stereo) ... dump of sz links" << std::endl;
350 Dump( forLine, "sort flag", " " );
351#endif
352
353 //...Check # of sz links...
354 if ( forLine.length() < _lineSelector.nLinksStereo() )
355 {
356#ifdef TRKRECO_DEBUG_DETAIL
357 std::cout << _name << "(stereo) ... rejected by sz nLinks(";
358 std::cout << forLine.length() << ") < ";
359 std::cout << _lineSelector.nLinks() << std::endl;
360#endif
361 HepAListDeleteAll( forLine );
362 return NULL;
363 }
364
365 //...Make a line...
366 unsigned nLine = forLine.length();
367 TLine0 line( forLine );
368 int err = line.fit2sp();
369 if ( err < 0 )
370 {
371 err = line.fit2p();
372 if ( err < 0 ) err = line.fit();
373 }
374
375 //...Linear fit...
376 if ( err < 0 )
377 {
378#ifdef TRKRECO_DEBUG_DETAIL
379 std::cout << _name << "(stereo) ... linear fit failure. nLinks(";
380 std::cout << forLine.length() << ")" << std::endl;
381#endif
382 HepAListDeleteAll( forLine );
383 return NULL;
384 }
385
386#ifdef TRKRECO_DEBUG_DETAIL
387 std::cout << _name << "(stereo) ... dump of left-right" << std::endl;
388#endif
389
390 //...Decide Left or Right...
391 AList<TMLink> forNewLine;
392 for ( unsigned i = 0; i < nLine; i++ )
393 {
394 TMLink* t = forLine[i];
395 TMLink* tl = new TMLink( *t );
396 TMLink* tr = new TMLink( *t );
397
398 tl->leftRight( WireHitLeft );
399 tr->leftRight( WireHitRight );
400
401 int err = track.szPosition( *tl );
402 if ( err )
403 {
404 delete tl;
405 tl = NULL;
406 }
407 err = track.szPosition( *tr );
408 if ( err )
409 {
410 delete tr;
411 tr = NULL;
412 }
413 if ( ( tl == NULL ) && ( tr == NULL ) ) continue;
414
415 TMLink* best;
416 if ( tl == NULL ) best = tr;
417 else if ( tr == NULL ) best = tl;
418 else
419 {
420 if ( line.distance( *tl ) < line.distance( *tr ) )
421 {
422 best = tl;
423 delete tr;
424 }
425 else
426 {
427 best = tr;
428 delete tl;
429 }
430 }
431
432#ifdef TRKRECO_DEBUG_DETAIL
433 std::cout << " ";
434 std::cout << t->wire()->layerId() << "-";
435 std::cout << t->wire()->localId();
436 if ( tl != NULL ) std::cout << ",left " << tl->position() << "," << line.distance( *tl );
437 if ( tr != NULL ) std::cout << ",right " << tr->position() << "," << line.distance( *tr );
438 std::cout << std::endl;
439#endif
440
441 best->link( t->link() );
442 forNewLine.append( best );
443 }
444
445 //...Check # of sz links...
446 if ( forNewLine.length() < _lineSelector.nLinksStereo() )
447 {
448#ifdef TRKRECO_DEBUG_DETAIL
449 std::cout << _name << "(stereo) ... rejected by lr nLinks(";
450 std::cout << forNewLine.length() << ") < ";
451 std::cout << _lineSelector.nLinks() << std::endl;
452#endif
453 HepAListDeleteAll( forLine );
454 HepAListDeleteAll( forNewLine );
455 return NULL;
456 }
457
458 //...Create new line...
459#ifdef TRKRECO_DEBUG_DETAIL
460 std::cout << _name << "(stereo) ... creating a new line" << std::endl;
461#endif
462 unsigned nNewLine = forNewLine.length();
463 TLine0 newLine( forNewLine );
464
465 //... Remove extremely bad points
466 newLine.removeChits();
467
468 //...Make a seed track again
469 err = newLine.fit2sp();
470 if ( err < 0 )
471 {
472 err = newLine.fit2p();
473 if ( err < 0 ) err = newLine.fit();
474 }
475
476 //...Linear fit...
477 if ( err < 0 )
478 {
479#ifdef TRKRECO_DEBUG_DETAIL
480 std::cout << _name << "(stereo) ... 2nd linear fit failure. nLinks(";
481 std::cout << forNewLine.length() << ")" << std::endl;
482#endif
483 HepAListDeleteAll( forLine );
484 HepAListDeleteAll( forNewLine );
485 return NULL;
486 }
487
488 //...Remove bad points...
489 AList<TMLink> bad;
490 newLine.refine( bad, 40. );
491 err = newLine.fit();
492 newLine.refine( bad, 30. );
493 err = newLine.fit();
494 newLine.refine( bad, 20. );
495 err = newLine.fit();
496 newLine.refine( bad, 10. );
497 err = newLine.fit();
498 float R = fabs( track.helix().curv() );
499 if ( R > 80. )
500 {
501 newLine.refine( bad, 5. );
502 err = newLine.fit();
503 }
504
505 //...Linear fit again...
506 if ( err < 0 )
507 {
508 HepAListDeleteAll( forLine );
509 HepAListDeleteAll( forNewLine );
510#ifdef TRKRECO_DEBUG_DETAIL
511 std::cout << " appendStereo cut ... new line 2nd linear fit failure. ";
512 std::cout << "# of links = " << n << "," << nLine;
513 std::cout << "," << nNewLine << std::endl;
514#endif
515 return NULL;
516 }
517
518 //...3D fit...
519 const AList<TMLink>& good = newLine.links();
520 unsigned nn = good.length();
521 for ( unsigned i = 0; i < nn; i++ ) { track.append( *good[i]->link() ); }
522 Vector a( 5 );
523 a = track.helix().a();
524 a[4] = track.charge() * newLine.a();
525 track._helix->a( a );
526
527 //...Refine...
528 err = fit( track );
529 track.refine( bad, _trackSelector.maxSigma() * 100. );
530 err = fit( track );
531 track.refine( bad, _trackSelector.maxSigma() * 10. );
532 err = fit( track );
533 track.refine( bad, _trackSelector.maxSigma() );
534 err = fit( track );
535
536 //...Test it...
537 if ( !_trackSelector.select( track ) )
538 {
539 HepAListDeleteAll( forLine );
540 HepAListDeleteAll( forNewLine );
541 return NULL;
542 }
543
544 //...Termination...
545 HepAListDeleteAll( forLine );
546 HepAListDeleteAll( forNewLine );
547 return &track;
548}
549
550TTrack* TBuilder0::buildStereo( TTrack& track, const AList<TMLink>& list ) const {
551#ifdef TRKRECO_DEBUG_DETAIL
552 std::cout << _name << "(stereo) ... building in 3D" << std::endl;
553 std::cout << "... dump of stereo candidate hits" << std::endl;
554 Dump( list, "sort flag", " " );
555#endif
556
557 //...Check # of links...
558 if ( list.length() < _lineSelector.nLinksStereo() )
559 {
560#ifdef TRKRECO_DEBUG_DETAIL
561 std::cout << "... rejected by nLinks(";
562 std::cout << list.length() << ") < ";
563 std::cout << _lineSelector.nLinks() << std::endl;
564#endif
565 return NULL;
566 }
567
568 //...Setup...
569 int ichg;
570 if ( track.helix().curv() > 0. ) ichg = 1;
571 else ichg = -1;
572 unsigned nlyr[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
573 unsigned llyr[18] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
574
575 //...Check
576 //... # of hits in a wire,
577 //... 2 contineus hits or not
578 //...
579 //...and Calculate s and z for every links...
580 unsigned n = list.length();
581 AList<TMLink> forLine;
582 unsigned iforLine = 0;
583 for ( unsigned i = 0; i < n; i++ )
584 {
585 TMLink* l = list[i];
586
587 if ( i < n - 1 )
588 {
589 TMLink* lnext = list[i + 1];
590
591 //...2 consective hits ?...
592 // int LorR = check2CnHits(* l, * lnext, ichg);
593 int LorR = consectiveHits( *l, *lnext, ichg );
594 //...For debug...
595 // if (check2CnHits(* l, * lnext, ichg) !=
596 // consectiveHits(* l, * lnext, ichg)) {
597 // AList<TMLink> tmp;
598 // tmp.append(l);
599 // tmp.append(lnext);
600 // std::cout << "!!! consective diff !!!" << std::endl;
601 // std::cout << " original = " << check2CnHits(* l, * lnext, ichg);
602 // std::cout << ", iw = " << consectiveHits(* l, * lnext, ichg);
603 // std::cout << std::endl;
604 // Dump(tmp, "detail");
605 // }
606 //...For debug end...
607
608 //...Left/right solved by two consective hits...
609 if ( LorR == WireHitLeft || LorR == WireHitRight )
610 {
611
612 //... Set Left/Right
613 if ( LorR == WireHitLeft )
614 {
615 l->leftRight( WireHitLeft );
616 lnext->leftRight( WireHitRight );
617 }
618 else
619 {
620 l->leftRight( WireHitRight );
621 lnext->leftRight( WireHitLeft );
622 }
623
624 //...Calculate z...
625 int err1 = track.szPosition( *l );
626 int err2 = track.szPosition( *lnext );
627 if ( err1 == 0 && err2 == 0 )
628 {
629 double deltaZ = fabs( l->position().y() - lnext->position().y() );
630 if ( deltaZ < 1.5 )
631 {
632
633 //... O.K. l and lnext should be good 2 consectvie hits
634 l->zStatus( 20 );
635 lnext->zStatus( 20 );
636 l->zPair( iforLine + 1 );
637 lnext->zPair( iforLine );
638 }
639#ifdef TRKRECO_DEBUG_DETAIL
640 else
641 {
642 std::cout << " ... rejected because delta z > 1.5";
643 std::cout << std::endl;
644 }
645#endif
646 }
647#ifdef TRKRECO_DEBUG_DETAIL
648 else
649 {
650 if ( err1 )
651 {
652 std::cout << " ... s-z calculation error with ";
653 std::cout << l->wire()->name() << std::endl;
654 }
655 if ( err2 )
656 {
657 std::cout << " ... s-z calculation error with ";
658 std::cout << lnext->wire()->name() << std::endl;
659 }
660 }
661#endif
662 }
663 }
664
665 //... Calculate s and z...
666 //... Aleady solved
667 if ( l->zStatus() == 20 )
668 {
669 TMLink* t = new TMLink( *l );
670 t->zStatus( l->zStatus() );
671 t->zPair( l->zPair() );
672 int err = track.szPosition( *t );
673 if ( err )
674 {
675 delete t;
676 continue;
677 }
678 else
679 {
680 //...Store the sz link...
681 t->link( l );
682 forLine.append( t );
683 //... Check # of hits in a wire layer and Clustering them
684 // nlyr[layertoStlayer(l->wire()->layerId())] += 1;
685 nlyr[l->wire()->layer()->axialStereoLayerId()] += 1;
686 // llyr[layertoStlayer(l->wire()->layerId())] = iforLine;
687 llyr[l->wire()->layer()->axialStereoLayerId()] = iforLine;
688 iforLine += 1;
689 }
690 }
691 else
692 {
693 //...Assuming wire position...
694 //... for left
695 TMLink* tl = new TMLink( *l );
696 tl->leftRight( WireHitLeft );
697 tl->zStatus( -10 );
698 tl->zPair( 0 );
699 int err = track.szPosition( *tl );
700 if ( err ) { delete tl; }
701 else
702 {
703 //...Store the sz link...
704 tl->link( l );
705 forLine.append( tl );
706 //... Check # of hits in a wire layer and Clustering them
707 // nlyr[layertoStlayer(l->wire()->layerId())] += 1;
708 nlyr[l->wire()->layer()->axialStereoLayerId()] += 1;
709 // llyr[layertoStlayer(l->wire()->layerId())] = iforLine;
710 llyr[l->wire()->layer()->axialStereoLayerId()] = iforLine;
711 //...
712 iforLine += 1;
713 }
714
715 //... for right
716 TMLink* tr = new TMLink( *l );
717 tr->leftRight( WireHitRight );
718 tr->zStatus( -10 );
719 tr->zPair( 0 );
720 err = track.szPosition( *tr );
721 if ( err ) { delete tr; }
722 else
723 {
724 //...Store the sz link...
725 tr->link( l );
726 forLine.append( tr );
727 //... Check # of hits in a wire layer and Clustering them
728 // nlyr[layertoStlayer(l->wire()->layerId())] += 1;
729 nlyr[l->wire()->layer()->axialStereoLayerId()] += 1;
730 // llyr[layertoStlayer(l->wire()->layerId())] = iforLine;
731 llyr[l->wire()->layer()->axialStereoLayerId()] = iforLine;
732 iforLine += 1;
733 }
734 }
735 }
736
737#ifdef TRKRECO_DEBUG_DETAIL
738 std::cout << "... dump of sz links" << std::endl;
739 Dump( forLine, "sort flag", " " );
740 for ( unsigned i = 0; i < 18; i++ )
741 { std::cout << i << " : " << nlyr[i] << ", " << llyr[i] << std::endl; }
742#endif
743
744 //...Check # of sz links...
745 if ( forLine.length() < _lineSelector.nLinksStereo() )
746 {
747#ifdef TRKRECO_DEBUG_DETAIL
748 std::cout << "... rejected by sz nLinks(";
749 std::cout << forLine.length() << ") < ";
750 std::cout << _lineSelector.nLinks() << std::endl;
751#endif
752 HepAListDeleteAll( forLine );
753 return NULL;
754 }
755
756 //... Select the best segment in each Superlayer
757 AList<TMLink> gdSLine0;
758 AList<TMLink> gdSLine1;
759 AList<TMLink> gdSLine2;
760 AList<TMLink> gdSLine3;
761 AList<TMLink> gdSLine4;
762 double min_chi2[5] = { 9999., 9999., 9999., 9999., 9999. };
763 double min_a[5] = { 9999., 9999., 9999., 9999., 9999. };
764
765#ifdef TRKRECO_DEBUG_DETAIL
766 // TTrackBase base;
767 // base.append(forLine);
768 // TWindow baseWindow("s-z window");
769 // baseWindow.append(base);
770
771 bool display = false;
772// for (unsigned i = 0; i < forLine.length(); i++) {
773// if (forLine[i]->wire()->name() == "15=31")
774// display = true;
775// }
776#endif
777
778 //...For all stereo super layers...
779 for ( unsigned isl = 0; isl < 5; isl++ )
780 {
781 AList<TMLink> tmpLine;
782 AList<TMLink> goodLine;
783
784#ifdef TRKRECO_DEBUG_DETAIL
785 std::cout << " ... stereo supter layer " << isl << std::endl;
786#endif
787
788 //...3-stereo-layer case...
789 if ( isl < 2 )
790 {
791
792 //...Initialize...
793 unsigned ily1 =
794 TMDC::getTMDC( "1.0" )->superLayer( isl * 2 + 1 )->first()->axialStereoLayerId();
795 unsigned ily2 = ily1 + 1;
796 unsigned ily3 = ily2 + 1;
797
798 //...Loop for the first layer...
799 unsigned ilst = llyr[ily1] + 1;
800 unsigned ifst = ilst - nlyr[ily1];
801 for ( unsigned i = ifst; i < ilst; i++ )
802 {
803 TMLink& l = *forLine[i];
804 tmpLine.append( l );
805 TMLink& m = *forLine[l.zPair()];
806 if ( l.zStatus() == 20 )
807 {
808 if ( l.zPair() < i ) { tmpLine.append( m ); }
809 else
810 {
811 //...Already considered...
812 tmpLine.remove( l );
813 continue;
814 }
815 }
816
817 //...Loop for the second layer...
818 unsigned jlst = llyr[ily2] + 1;
819 unsigned jfst = jlst - nlyr[ily2];
820 for ( unsigned j = jfst; j < jlst; j++ )
821 {
822 TMLink& l2 = *forLine[j];
823 tmpLine.append( l2 );
824 TMLink& m2 = *forLine[l2.zPair()];
825 if ( l2.zStatus() == 20 )
826 {
827 if ( l2.zPair() < j ) { tmpLine.append( m2 ); }
828 else
829 {
830 //...Already considered.
831 tmpLine.remove( l2 );
832 continue;
833 }
834 }
835
836 //...Loop for the third layer...
837 unsigned klst = llyr[ily3] + 1;
838 unsigned kfst = klst - nlyr[ily3];
839 for ( unsigned k = kfst; k < klst; k++ )
840 {
841 TMLink& l3 = *forLine[k];
842 tmpLine.append( l3 );
843 TMLink& m3 = *forLine[l3.zPair()];
844 if ( l3.zStatus() == 20 )
845 {
846 if ( l3.zPair() < k ) { tmpLine.append( m3 ); }
847 else
848 {
849 //...Already considered.
850 tmpLine.remove( l3 );
851 continue;
852 }
853 }
854
855 //... Check the hits in neighbor wirelayer
856 //...
857 //... x|o|o|x |
858 //... |o| |
859 //... V IP , o means OK.
860 //...
861 int relation12 = -1;
862 int relation23 = -1;
863
864 //... Check the hit in fist and it in the second layer
865 relation12 =
866 checkHits( l.hit()->wire()->localId(), l2.hit()->wire()->localId(), isl );
867 //...
868 if ( l.zStatus() == 20 && relation12 < 0 )
869 relation12 =
870 checkHits( m.hit()->wire()->localId(), l2.hit()->wire()->localId(), isl );
871
872 //...Check the hit in second and it in the third layer
873 relation23 =
874 checkHits( l2.hit()->wire()->localId(), l3.hit()->wire()->localId(), isl );
875 if ( l.zStatus() == 20 && relation23 < 0 )
876 relation23 =
877 checkHits( m2.hit()->wire()->localId(), l3.hit()->wire()->localId(), isl );
878
879 //...Bad relation...
880 if ( relation12 || relation23 )
881 {
882#ifdef TRKRECO_DEBUG_DETAIL
883 std::cout << " ... bad relations";
884 std::cout << " : segment rejected";
885 Dump( tmpLine, "detail stereo", " " );
886 if ( display )
887 {
888 TLine0 line( tmpLine );
889 int err = line.fit();
890 // baseWindow.append(line);
891 // baseWindow.draw();
892 // baseWindow.wait();
893 // baseWindow.remove(line);
894 }
895#endif
896
897 tmpLine.remove( l3 );
898 if ( l3.zStatus() == 20 ) tmpLine.remove( m3 );
899 continue;
900 }
901
902 //...Make a segment
903 unsigned ntmp = tmpLine.length();
904 if ( ntmp < 3 )
905 { //...Is this needed???...
906 std::cout << "!!! is this possible !!!???" << std::endl;
907
908 tmpLine.remove( l3 );
909 if ( l3.zStatus() == 20 ) tmpLine.remove( m3 );
910 continue;
911 }
912
913 //...Do the linear fit
914 TLine0 line( tmpLine );
915 int err = line.fit();
916 if ( err < 0 )
917 {
918#ifdef TRKRECO_DEBUG_DETAIL
919 std::cout << " ... line fit error";
920 std::cout << " : segment rejected" << std::endl;
921 Dump( line.links(), "detail stereo", " " );
922 if ( display )
923 {
924 // baseWindow.append(line);
925 // baseWindow.draw();
926 // baseWindow.wait();
927 // baseWindow.remove(line);
928 }
929#endif
930
931 tmpLine.remove( l3 );
932 if ( l3.zStatus() == 20 ) tmpLine.remove( m3 );
933 continue;
934 }
935
936 //...Check quality for the fit
937 //... if |z0| > 70 or chi2 > 0.7, remove it
938 // double chi2 = line.chi2()/(ntmp-2);
939 double chi2 = line.reducedChi2();
940 if ( fabs( line.b() ) < _stereoZ3 && chi2 < _stereoChisq3 && chi2 < min_chi2[isl] )
941 {
942 goodLine = tmpLine;
943 min_chi2[isl] = chi2;
944 min_a[isl] = line.a();
945#ifdef TRKRECO_DEBUG_DETAIL
946 std::cout << " ... segment accepted : " << std::endl;
947#endif
948 }
949#ifdef TRKRECO_DEBUG_DETAIL
950 else
951 {
952 std::cout << " ... bad quality : segment rejected :";
953 std::cout << "chi2=" << chi2;
954 std::cout << ", b=" << line.b() << std::endl;
955 }
956 Dump( line.links(), "detail stereo", " " );
957 if ( display )
958 {
959 // baseWindow.append(line);
960 // baseWindow.draw();
961 // baseWindow.wait();
962 // baseWindow.remove(line);
963 }
964#endif
965
966 //...end of third loop
967 if ( l3.zStatus() == 20 ) tmpLine.remove( m3 );
968 tmpLine.remove( l3 );
969 }
970 //... end of second loop
971 if ( l2.zStatus() == 20 ) tmpLine.remove( m2 );
972 tmpLine.remove( l2 );
973 }
974 //... end of first loop
975 if ( l.zStatus() == 20 ) tmpLine.remove( m );
976 tmpLine.remove( l );
977 }
978
979 //... Keep best segments
980 if ( isl == 0 ) gdSLine0 = goodLine;
981 if ( isl == 1 ) gdSLine1 = goodLine;
982 goodLine.removeAll();
983 }
984 else
985 {
986
987 //... Nlayer == 4
988
989 //...Initialize
990 // unsigned ily1 = firstStlayer(isl);
991 unsigned ily1 =
992 TMDC::getTMDC( "1.0" )->superLayer( isl * 2 + 1 )->first()->axialStereoLayerId();
993 unsigned ily2 = ily1 + 1;
994 unsigned ily3 = ily2 + 1;
995 unsigned ily4 = ily3 + 1;
996 if ( nlyr[ily1] == 0 || nlyr[ily2] == 0 || nlyr[ily3] == 0 || nlyr[ily4] == 0 ) continue;
997
998 //...loop for the first layer
999 unsigned ilst = llyr[ily1] + 1;
1000 unsigned ifst = ilst - nlyr[ily1];
1001 for ( unsigned i = ifst; i < ilst; i++ )
1002 {
1003 TMLink& l = *forLine[i];
1004 tmpLine.append( l );
1005 TMLink& m = *forLine[l.zPair()];
1006 if ( l.zStatus() == 20 )
1007 {
1008 if ( l.zPair() < i ) { tmpLine.append( m ); }
1009 else
1010 {
1011 tmpLine.remove( l );
1012 continue;
1013 }
1014 }
1015
1016 //...loop for the second layer
1017 unsigned jlst = llyr[ily2] + 1;
1018 unsigned jfst = jlst - nlyr[ily2];
1019 for ( unsigned j = jfst; j < jlst; j++ )
1020 {
1021 TMLink& l2 = *forLine[j];
1022 tmpLine.append( l2 );
1023 TMLink& m2 = *forLine[l2.zPair()];
1024 if ( l2.zStatus() == 20 )
1025 {
1026 if ( l2.zPair() < j ) { tmpLine.append( m2 ); }
1027 else
1028 {
1029 tmpLine.remove( l2 );
1030 continue;
1031 }
1032 }
1033
1034 //...loop for the third layer
1035 unsigned klst = llyr[ily3] + 1;
1036 unsigned kfst = klst - nlyr[ily3];
1037 for ( unsigned k = kfst; k < klst; k++ )
1038 {
1039 TMLink& l3 = *forLine[k];
1040 tmpLine.append( l3 );
1041 TMLink& m3 = *forLine[l3.zPair()];
1042 if ( l3.zStatus() == 20 )
1043 {
1044 if ( l3.zPair() < k ) { tmpLine.append( m3 ); }
1045 else
1046 {
1047 tmpLine.remove( l3 );
1048 continue;
1049 }
1050 }
1051
1052 //...loop for the 4th layer
1053 unsigned hlst = llyr[ily4] + 1;
1054 unsigned hfst = hlst - nlyr[ily4];
1055 for ( unsigned h = hfst; h < hlst; h++ )
1056 {
1057 TMLink& l4 = *forLine[h];
1058 tmpLine.append( l4 );
1059 TMLink& m4 = *forLine[l4.zPair()];
1060 if ( l4.zStatus() == 20 )
1061 {
1062 if ( l4.zPair() < h ) { tmpLine.append( m4 ); }
1063 else
1064 {
1065 tmpLine.remove( l4 );
1066 continue;
1067 }
1068 }
1069
1070 //...Check the relation between the hit
1071 // in neighbor wirelayer
1072 //...
1073 //... x|o|o|x |
1074 //... |o| |
1075 //... V IP , o means OK.
1076 //...
1077 int relation12 = -1;
1078 int relation23 = -1;
1079 int relation34 = -1;
1080
1081 //...For debug...
1082 // if (l.hit()->wire()->consective(* l2.hit()->wire())
1083 //...For debug end...
1084
1085 //... Check the hit in fist and it in the second
1086 relation12 =
1087 checkHits( l.hit()->wire()->localId(), l2.hit()->wire()->localId(), isl );
1088 if ( l.zStatus() == 20 && relation12 < 0 )
1089 relation12 =
1090 checkHits( m.hit()->wire()->localId(), l2.hit()->wire()->localId(), isl );
1091
1092 //... Check the hit in second and it in the third
1093 relation23 =
1094 checkHits( l2.hit()->wire()->localId(), l3.hit()->wire()->localId(), isl );
1095 if ( l.zStatus() == 20 && relation23 < 0 )
1096 relation23 =
1097 checkHits( m2.hit()->wire()->localId(), l3.hit()->wire()->localId(), isl );
1098
1099 //... Check the hit in second and it in the forth
1100 relation34 =
1101 checkHits( l3.hit()->wire()->localId(), l4.hit()->wire()->localId(), isl );
1102 if ( l3.zStatus() == 20 && relation34 < 0 )
1103 relation34 =
1104 checkHits( m3.hit()->wire()->localId(), l4.hit()->wire()->localId(), isl );
1105
1106 //...remove Bad segments
1107 if ( relation12 || relation23 || relation34 )
1108 {
1109
1110#ifdef TRKRECO_DEBUG_DETAIL
1111 std::cout << " ... bad relations";
1112 std::cout << " : segment rejected";
1113 Dump( tmpLine, "detail stereo", " " );
1114 if ( display )
1115 {
1116 TLine0 line( tmpLine );
1117 int err = line.fit();
1118 // baseWindow.append(line);
1119 // baseWindow.draw();
1120 // baseWindow.wait();
1121 // baseWindow.remove(line);
1122 }
1123#endif
1124
1125 tmpLine.remove( l4 );
1126 if ( l4.zStatus() == 20 ) tmpLine.remove( m4 );
1127 continue;
1128 }
1129
1130 //...Make a segment
1131 unsigned ntmp = tmpLine.length();
1132 if ( ntmp < 4 )
1133 {
1134 tmpLine.remove( l4 );
1135 if ( l4.zStatus() == 20 ) tmpLine.remove( m4 );
1136 continue;
1137 }
1138
1139 //...Do the linear fit
1140 TLine0 line( tmpLine );
1141 int err = line.fit();
1142 if ( err < 0 )
1143 {
1144#ifdef TRKRECO_DEBUG_DETAIL
1145 std::cout << " ... line fit error";
1146 std::cout << " : segment rejected" << std::endl;
1147 Dump( line.links(), "detail stereo", " " );
1148 if ( display )
1149 {
1150 // baseWindow.append(line);
1151 // baseWindow.draw();
1152 // baseWindow.wait();
1153 // baseWindow.remove(line);
1154 }
1155#endif
1156
1157 tmpLine.remove( l4 );
1158 if ( l4.zStatus() == 20 ) tmpLine.remove( m4 );
1159 continue;
1160 }
1161
1162 //...Check qualit of the fit
1163 //... if |z0| > 20 or chi2 > 0.5, remove it
1164 double chi2 = line.reducedChi2();
1165 if ( fabs( line.b() ) < _stereoZ4 && chi2 < _stereoChisq4 &&
1166 chi2 < min_chi2[isl] )
1167 {
1168
1169 //...Keep good segments
1170 goodLine = tmpLine;
1171 min_chi2[isl] = chi2;
1172 min_a[isl] = line.a();
1173
1174#ifdef TRKRECO_DEBUG_DETAIL
1175 std::cout << " segment accepted : " << std::endl;
1176#endif
1177 }
1178#ifdef TRKRECO_DEBUG_DETAIL
1179 else
1180 {
1181 std::cout << " ... bad quality : segment rejected :";
1182 std::cout << " chi2=" << chi2;
1183 std::cout << ", b=" << line.b() << std::endl;
1184 }
1185 Dump( line.links(), "detail stereo", " " );
1186 if ( display )
1187 {
1188 // baseWindow.append(line);
1189 // baseWindow.draw();
1190 // baseWindow.wait();
1191 // baseWindow.remove(line);
1192 }
1193#endif
1194
1195 //...end of the 4th loop
1196 if ( l4.zStatus() == 20 ) tmpLine.remove( m4 );
1197 tmpLine.remove( l4 );
1198 }
1199 //...end of the third loop
1200 if ( l3.zStatus() == 20 ) tmpLine.remove( m3 );
1201 tmpLine.remove( l3 );
1202 }
1203 //...end of the second loop
1204 if ( l2.zStatus() == 20 ) tmpLine.remove( m2 );
1205 tmpLine.remove( l2 );
1206 }
1207 //... end of the first loop
1208 if ( l.zStatus() == 20 ) tmpLine.remove( m );
1209 tmpLine.remove( l );
1210 }
1211
1212 //...
1213 if ( isl == 2 ) gdSLine2 = goodLine;
1214 if ( isl == 3 ) gdSLine3 = goodLine;
1215 if ( isl == 4 ) gdSLine4 = goodLine;
1216 goodLine.removeAll();
1217 }
1218 }
1219
1220 //... Link segments
1221
1222 //...Check how many segments are made
1223 unsigned Nsgmnts[5] = { 0, 0, 0, 0, 0 };
1224 Nsgmnts[0] = gdSLine0.length();
1225 Nsgmnts[1] = gdSLine1.length();
1226 Nsgmnts[2] = gdSLine2.length();
1227 Nsgmnts[3] = gdSLine3.length();
1228 Nsgmnts[4] = gdSLine4.length();
1229
1230 unsigned NusedSgmnts = 0;
1231 for ( unsigned jsl = 0; jsl < 5; jsl++ )
1232 {
1233 if ( Nsgmnts[jsl] > 0 ) NusedSgmnts += 1;
1234 }
1235
1236 //...Require at least one Segment
1237 if ( NusedSgmnts == 0 )
1238 {
1239#ifdef TRKRECO_DEBUG_DETAIL
1240 std::cout << "... rejected because no segment found" << std::endl;
1241#endif
1242 HepAListDeleteAll( forLine );
1243 return NULL;
1244 }
1245
1246 //... Make a line
1247 AList<TMLink> forNewLine;
1248 forNewLine.append( gdSLine0 );
1249 forNewLine.append( gdSLine1 );
1250 forNewLine.append( gdSLine2 );
1251 forNewLine.append( gdSLine3 );
1252 forNewLine.append( gdSLine4 );
1253
1254 //...
1255 unsigned nNewLine = forNewLine.length();
1256 float R = fabs( track.helix().curv() );
1257 TLine0 newLine( forNewLine );
1258
1259 //...Do linear fit
1260 int err = newLine.fit();
1261 // double newLine_chi2 = newLine.chi2()/(nNewLine - 2);
1262 double newLine_chi2 = newLine.reducedChi2();
1263 double newLine_a = newLine.a();
1264
1265 //...Check the quality of the line.
1266
1267 //...If chi2 > 0.25 for R > 80.( > 0.8 for R < 80.), refine the line.
1268 // if(((R > 80. && newLine_chi2 > 0.25) ||(R < 80. && newLine_chi2 > 0.8))&& NusedSgmnts >
1269 // 1){
1270 if ( ( ( R > 80. && newLine_chi2 > 4.0 ) || ( R < 80. && newLine_chi2 > 13.0 ) ) &&
1271 NusedSgmnts > 1 )
1272 {
1273 //...Look at difference between the slope of the line and that of segment
1274 double max_diff_a = 0.;
1275 unsigned this_sly = 999;
1276 for ( unsigned isl = 0; isl < 5; isl++ )
1277 {
1278 if ( Nsgmnts[isl] == 0 ) continue;
1279 double diff_a = fabs( ( min_a[isl] - newLine_a ) / newLine_a );
1280 if ( diff_a > max_diff_a )
1281 {
1282 max_diff_a = diff_a;
1283 this_sly = isl;
1284 }
1285 }
1286
1287 //...If max slope diff. > 0.4 for R < 50.(> 0.3 for R < 50), remove it.
1288 if ( ( R < 50. && max_diff_a > 0.4 ) || ( R > 50. && max_diff_a > 0.3 ) )
1289 {
1290
1291 //...clear # of entries in the segement
1292 Nsgmnts[this_sly] = 0;
1293
1294 //... remove the worst setment
1295 if ( this_sly == 0 ) { newLine.removeSLY( gdSLine0 ); }
1296 else if ( this_sly == 1 ) { newLine.removeSLY( gdSLine1 ); }
1297 else if ( this_sly == 2 ) { newLine.removeSLY( gdSLine2 ); }
1298 else if ( this_sly == 3 ) { newLine.removeSLY( gdSLine3 ); }
1299 else if ( this_sly == 4 ) { newLine.removeSLY( gdSLine4 ); }
1300
1301 //... fit again
1302 unsigned NnewLine_2 = newLine.links().length();
1303 if ( NnewLine_2 < 3 )
1304 {
1305#ifdef TRKRECO_DEBUG_DETAIL
1306 std::cout << "... rejected because of few links for line" << std::endl;
1307#endif
1308 HepAListDeleteAll( forLine );
1309 return NULL;
1310 }
1311
1312 int err = newLine.fit();
1313 if ( err < 0 )
1314 {
1315#ifdef TRKRECO_DEBUG_DETAIL
1316 std::cout << "... rejected because of line fit failure" << std::endl;
1317#endif
1318 HepAListDeleteAll( forLine );
1319 return NULL;
1320 }
1321 double newLine_chi2_2 = newLine.chi2() / NnewLine_2;
1322 }
1323 }
1324
1325 //...Remove bad points...
1326 AList<TMLink> bad;
1327 double maxSigma = 1.0;
1328 if ( R < 80 ) maxSigma = 1.5;
1329 newLine.refine( bad, maxSigma );
1330 if ( newLine.links().length() < 2 )
1331 {
1332#ifdef TRKRECO_DEBUG_DETAIL
1333 std::cout << "... rejected because of few links for line after refine";
1334 std::cout << std::endl;
1335#endif
1336 HepAListDeleteAll( forLine );
1337 return NULL;
1338 }
1339 err = newLine.fit();
1340 if ( err < 0 )
1341 {
1342#ifdef TRKRECO_DEBUG_DETAIL
1343 std::cout << "... rejected because of line fit failure(2)" << std::endl;
1344#endif
1345 HepAListDeleteAll( forLine );
1346 return NULL;
1347 }
1348
1349 //...Append hits, if the distance between the line and hits <
1350 for ( unsigned isl = 0; isl < 5; isl++ )
1351 {
1352 if ( Nsgmnts[isl] == 0 )
1353 {
1354 double maxdist = 0.5;
1355 if ( R < 80 ) maxdist = 1.25;
1356 newLine.appendByszdistance( forLine, 2 * isl + 1, maxdist );
1357 }
1358 }
1359 err = newLine.fit();
1360 if ( err < 0 )
1361 {
1362#ifdef TRKRECO_DEBUG_DETAIL
1363 std::cout << "... rejected because of line fit failure(3)" << std::endl;
1364#endif
1365 HepAListDeleteAll( forLine );
1366 return NULL;
1367 }
1368
1369 //...Remove bad points again...
1370 maxSigma = 1.0;
1371 newLine.refine( bad, maxSigma );
1372 if ( newLine.links().length() < 2 )
1373 {
1374#ifdef TRKRECO_DEBUG_DETAIL
1375 std::cout << "... rejected because of few links for line after 2nd refine";
1376 std::cout << std::endl;
1377#endif
1378 HepAListDeleteAll( forLine );
1379 return NULL;
1380 }
1381
1382 //...Linear fit again...
1383 err = newLine.fit();
1384 if ( err < 0 )
1385 {
1386 HepAListDeleteAll( forLine );
1387#ifdef TRKRECO_DEBUG_DETAIL
1388 std::cout << " appendStereo cut ... new line 2nd linear fit failure. ";
1389 std::cout << "# of links = " << n << ",";
1390 std::cout << "," << nNewLine << std::endl;
1391#endif
1392 return NULL;
1393 }
1394
1395 //...
1396 unsigned NnewLine_3 = newLine.links().length();
1397 double newLine_chi2_3 = newLine.chi2() / NnewLine_3;
1398
1399 //... Do we need quality cut?
1400 // if(newLine_chi2_3 > 10.) return NULL;
1401
1402 //...3D fit...
1403 const AList<TMLink>& good = newLine.links();
1404 unsigned nn = good.length();
1405 for ( unsigned i = 0; i < nn; i++ ) { track.append( *good[i]->link() ); }
1406
1407 //...Set initial values
1408 Vector a( 5 );
1409 a = track.helix().a();
1410 a[3] = newLine.b();
1411 a[4] = track.charge() * newLine.a();
1412 track._helix->a( a );
1413
1414 //...Refine...
1415 err = fit( track );
1416 if ( err < 0 ) goto discard;
1417 track.refine( bad, _stereoMaxSigma * 100. );
1418 err = fit( track );
1419 if ( err < 0 ) goto discard;
1420 track.refine( bad, _stereoMaxSigma * 10. );
1421 err = fit( track );
1422 if ( err < 0 ) goto discard;
1423 track.refine( bad, _stereoMaxSigma );
1424 err = fit( track );
1425 if ( err < 0 ) goto discard;
1426
1427 //...Test it...
1428 if ( !_trackSelector.select( track ) ) goto discard;
1429
1430 //...Termination...
1431 HepAListDeleteAll( forLine );
1432 return &track;
1433
1434 //...Something happened...
1435discard:
1436#ifdef TRKRECO_DEBUG_DETAIL
1437 std::cout << " ... rejected by fitting failure : ";
1438 std::cout << " err = " << err << std::endl;
1439#endif
1440 HepAListDeleteAll( forLine );
1441 return NULL;
1442}
1443
1444const TMSelector& TBuilder0::trackSelector( const TMSelector& a ) {
1445 _circleSelector.nLinks( a.nLinks() );
1446 _circleSelector.nSuperLayers( a.nSuperLayers() );
1447 _circleSelector.minPt( a.minPt() );
1448 _circleSelector.maxImpact( a.maxImpact() );
1449
1450 _trackSelector.nLinks( a.nLinks() );
1451 _trackSelector.nSuperLayers( a.nSuperLayers() );
1452 _trackSelector.maxSigma( a.maxSigma() );
1453
1454 _lineSelector.nLinksStereo( a.nLinksStereo() );
1455 _lineSelector.maxDistance( a.maxDistance() );
1456
1457 return a;
1458}
1459
1460void TBuilder0::appendClusters( TTrack& track, const AList<TMLink>& list ) const {
1461
1462 AList<TMLink> tmp = list;
1463
1464 //...Append them...
1465 track.appendByApproach( tmp, _trackSelector.maxSigma() * 2. / 3. );
1466
1467 //...Refine it...
1468 AList<TMLink> bad;
1469 fit( track );
1470 track.refine( bad, _trackSelector.maxSigma() );
1471}
1472
1473int TBuilder0::checkHits( unsigned i, unsigned j, unsigned isl ) const {
1474
1475 int nWr[5] = { 80, 128, 160, 208, 256 };
1476 int ilast = nWr[isl] - 1;
1477 int ilocal = (int)i;
1478 int jlocal = (int)j;
1479 //...
1480 if ( ilocal > 0 && ilocal < ilast )
1481 {
1482 if ( fabs( float( jlocal - ilocal ) ) > 1 ) { return -1; }
1483 else { return 0; }
1484 }
1485 else if ( ilocal == 0 )
1486 {
1487 if ( jlocal > 1 && jlocal < ilast ) { return -1; }
1488 else
1489 {
1490 if ( jlocal == ilast ) { return 0; }
1491 else if ( jlocal == 0 ) { return 0; }
1492 else if ( jlocal == 1 ) { return 0; }
1493 else { return -1; }
1494 }
1495 }
1496 else if ( ilocal == ilast )
1497 {
1498 if ( jlocal > 0 && jlocal < ilast - 1 ) { return -1; }
1499 else
1500 {
1501 if ( jlocal == ilast - 1 ) { return 0; }
1502 else if ( jlocal == ilast ) { return 0; }
1503 else if ( jlocal == 0 ) { return 0; }
1504 else { return -1; }
1505 }
1506 }
1507 //...
1508 return -1;
1509}
1510
1511void TBuilder0::salvage( TTrack& t, AList<TMLink>& hits ) const {
1512
1513 // if (t.type() == TrackTypeNormal) {
1514 salvageNormal( t, hits ); // Pt high enough; fromIP.
1515 return;
1516// }
1517// tmply for cosmic ray...
1518#ifdef TRKRECO_DEBUG_DETAIL
1519 std::cout << name() << " ... salvaging" << std::endl;
1520 std::cout << " # of given hits=" << hits.length() << std::endl;
1521#endif
1522
1523 unsigned nHits = hits.length();
1524 if ( nHits == 0 ) return;
1525
1526 //...Hit loop...
1527 AList<TMLink> candidates;
1528 for ( unsigned i = 0; i < nHits; i++ )
1529 {
1530 TMLink& l = *hits[i];
1531 const TMDCWireHit& h = *l.hit();
1532
1533 //...Already used?...
1534 if ( h.state() & WireHitUsed ) continue;
1535 candidates.append( l );
1536 }
1537
1538 //...Try to append this hit...
1539 // t.appendByApproach(candidates, 10.);
1540 t.appendByApproach( candidates, _salvageLevel );
1541 fit( t );
1542 hits.remove( candidates );
1543 t.assign( WireHitConformalFinder );
1544 t.finder( TrackOldConformalFinder );
1545 // t.assign(WireHitConformalFinder, TrackOldConformalFinder);
1546}
1547
1548void TBuilder0::salvageNormal( TTrack& t, AList<TMLink>& hits ) const {
1549
1550#ifdef TRKRECO_DEBUG_DETAIL
1551 std::cout << name() << " ... normal salvaging : ";
1552 std::cout << " # of hits=" << hits.length() << std::endl;
1553#endif
1554
1555 unsigned nHits = hits.length();
1556 if ( nHits == 0 ) return;
1557
1558 //...Determine direction to salvage...
1559 const HepPoint3D& center = t.helix().center();
1560 const HepVector3D a = ORIGIN - center;
1561
1562 //...Hit loop...
1563 AList<TMLink> candidates;
1564 for ( unsigned i = 0; i < nHits; i++ )
1565 {
1566 TMLink& l = *hits[i];
1567 const TMDCWireHit& h = *l.hit();
1568 if ( a.cross( h.xyPosition() ).z() * t.charge() > 0. ) continue;
1569
1570 //...Already used?...
1571 if ( h.state() & WireHitUsed ) continue;
1572 candidates.append( l );
1573 }
1574
1575 //...Try to append this hit...
1576 t.appendByApproach( candidates, 30. );
1577 fit( t );
1578 hits.remove( candidates );
1579 t.assign( WireHitConformalFinder );
1580 t.finder( TrackOldConformalFinder );
1581 // t.assign(WireHitConformalFinder, TrackOldConformalFinder);
1582}
1583
1584int TBuilder0::check2CnHits( TMLink& l, TMLink& s, int ichg ) const {
1585
1586 //...Check same layer ?...
1587 if ( l.hit()->wire()->layerId() != s.hit()->wire()->layerId() ) return -1;
1588
1589 //...Initialization...
1590 int nsl[11] = { 64, 80, 96, 128, 144, 160, 192, 208, 240, 256, 288 };
1591 float hcell[50] = { 0., 0., 0., 0., 0., 0., 0.687499, 0.747198, 0.806896,
1592 0., 0., 0., 0., 0., 0., 0.782967, 0.820598, 0.858229,
1593 0., 0., 0., 0., 0., 0.878423, 0.908646, 0.939845, 0.970068,
1594 0., 0., 0., 0., 0., 0.892908, 0.916188, 0.940219, 0.963499,
1595 0., 0., 0., 0., 0., 0.901912, 0.920841, 0.940382, 0.959312,
1596 0., 0., 0., 0., 0. };
1597
1598 int ilast = nsl[l.hit()->wire()->superLayerId()] - 1;
1599 int ilocal = l.hit()->wire()->localId();
1600 int jlocal = s.hit()->wire()->localId();
1601
1602 double ddist1 = l.hit()->drift() / hcell[l.hit()->wire()->layerId()];
1603 double ddist2 = s.hit()->drift() / hcell[s.hit()->wire()->layerId()];
1604 double ddist = 0.5 * ( ddist1 + ddist2 );
1605
1606 //...Case by case...
1607 if ( ilocal > 0 && ilocal < ilast )
1608 {
1609 if ( abs( jlocal - ilocal ) > 1 ) { return -20; }
1610 else
1611 {
1612 if ( ddist > 0.65 &&
1613 ( ( ddist1 > 0.7 && ddist2 > 0.7 ) || ( ddist1 > 0.95 || ddist2 > 0.95 ) ) )
1614 {
1615
1616 //...O.K. 2 consective hits
1617 if ( ichg > 0 ) { return 1; }
1618 else { return 0; }
1619 }
1620 else { return -20; }
1621 }
1622 }
1623 else if ( ilocal == 0 )
1624 {
1625 if ( jlocal > 1 && jlocal < ilast ) { return -20; }
1626 else
1627 {
1628 if ( ddist > 0.65 &&
1629 ( ( ddist1 > 0.7 && ddist2 > 0.7 ) || ( ddist1 > 0.95 || ddist2 > 0.95 ) ) )
1630 {
1631 if ( jlocal == ilast )
1632 {
1633
1634 //...O.K. 2 consective hits
1635 if ( ichg > 0 ) { return 0; }
1636 else { return 1; }
1637 }
1638 else if ( jlocal == 1 )
1639 {
1640 if ( ichg > 0 ) { return 1; }
1641 else { return 0; }
1642 }
1643 }
1644 }
1645 }
1646 else if ( ilocal == ilast )
1647 {
1648 if ( jlocal > 0 && jlocal < ilast - 1 ) { return -20; }
1649 else
1650 {
1651 if ( ddist > 0.65 &&
1652 ( ( ddist1 > 0.7 && ddist2 > 0.7 ) || ( ddist1 > 0.95 || ddist2 > 0.95 ) ) )
1653 {
1654 if ( jlocal == ilast - 1 )
1655 {
1656 //...O.K. 2 consective hits
1657 if ( ichg > 0 ) { return 0; }
1658 else { return 1; }
1659 }
1660 else if ( jlocal == 0 )
1661 {
1662 if ( ichg > 0 ) { return 1; }
1663 else { return 0; }
1664 }
1665 }
1666 else { return -20; }
1667 }
1668 }
1669
1670 //...fails
1671 return -50;
1672}
1673
1674int TBuilder0::consectiveHits( TMLink& l, TMLink& s, int ichg ) const {
1675
1676 //...Diagonal length of a cell...
1677 static float hcell[50] = { 0., 0., 0., 0., 0., 0., 0.687499, 0.747198, 0.806896,
1678 0., 0., 0., 0., 0., 0., 0.782967, 0.820598, 0.858229,
1679 0., 0., 0., 0., 0., 0.878423, 0.908646, 0.939845, 0.970068,
1680 0., 0., 0., 0., 0., 0.892908, 0.916188, 0.940219, 0.963499,
1681 0., 0., 0., 0., 0., 0.901912, 0.920841, 0.940382, 0.959312,
1682 0., 0., 0., 0., 0. };
1683 const TMDCWire& wire = *l.hit()->wire();
1684 const TMDCWire& next = *s.hit()->wire();
1685
1686 //...Check same layer ?...
1687 if ( wire.layerId() != next.layerId() ) return -1;
1688
1689 //...Are these consective ?...
1690 if ( !wire.consective( next ) ) return -20;
1691
1692 //...Drift distance...
1693 int ilast = wire.layer()->nWires() - 1;
1694 unsigned layerId = wire.layerId();
1695 double ddist1 = l.hit()->drift() / hcell[layerId];
1696 double ddist2 = s.hit()->drift() / hcell[layerId];
1697 double ddist = 0.5 * ( ddist1 + ddist2 );
1698
1699 //...Distance check...
1700 if ( ddist > 0.65 &&
1701 ( ( ddist1 > 0.7 && ddist2 > 0.7 ) || ( ddist1 > 0.95 || ddist2 > 0.95 ) ) )
1702 {
1703
1704 //...O.K. 2 consective hits
1705 if ( ichg > 0 ) return 1;
1706 else return 0;
1707 }
1708 else { return -50; }
1709}
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
tr
DOUBLE_PRECISION tr[3]
Definition EvtPhokharaDef.hh:63
s
XmlRpcServer s
Definition HelloServer.cpp:11
Vector
CLHEP::HepVector Vector
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:38
ORIGIN
const HepPoint3D ORIGIN
Constants.
Definition TMDCUtil.cxx:47
WireHitContinuous
#define WireHitContinuous
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:36
WireHitFittingValid
#define WireHitFittingValid
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:29
WireHitLeft
#define WireHitLeft
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:21
WireHitConformalFinder
#define WireHitConformalFinder
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:51
WireHitRight
#define WireHitRight
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:22
WireHitUsed
#define WireHitUsed
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:47
WireHitIsolated
#define WireHitIsolated
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:35
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
SameSuperLayer
AList< TMLink > SameSuperLayer(const AList< TMLink > &list, const TMLink &a)
returns links which are in the same super layer as 'a' or 'id'.
Definition TMLink.cxx:512
SuperLayer
unsigned SuperLayer(const AList< TMLink > &list)
returns super layer pattern.
Definition TMLink.cxx:565
NLayers
unsigned NLayers(const AList< TMLink > &links)
returns # of layers.
Definition TMLink.cxx:137
NSuperLayers
unsigned NSuperLayers(const AList< TMLink > &links)
returns # of layers.
Definition TMLink.cxx:582
TrackOldConformalFinder
#define TrackOldConformalFinder
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:23
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:217
TBuilder0::buildRphi
TTrack * buildRphi(const AList< TMLink > &) const
builds a r/phi track from TMLinks or from Segments.
Definition TBuilder0.cxx:72
TBuilder0::trackSelector
const TMSelector & trackSelector(void) const
returns a track selector.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TBuilder0.h:124
TBuilder0::~TBuilder0
virtual ~TBuilder0()
Destructor.
Definition TBuilder0.cxx:68
TBuilder0::buildStereo0
TTrack * buildStereo0(TTrack &track, const AList< TMLink > &) const
appends stereo hits to a track. (old version)
Definition TBuilder0.cxx:309
TBuilder0::dump
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TBuilder0.cxx:70
TBuilder0::name
const std::string & name(void) const
returns name.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TBuilder0.h:126
TBuilder0::TBuilder0
TBuilder0(const std::string &name)
Constructor.
Definition TBuilder0.cxx:31
TBuilder0::buildStereo
virtual TTrack * buildStereo(TTrack &track, const AList< TMLink > &) const
appends stereo hits to a track.
Definition TBuilder0.cxx:550
TBuilder0::salvage
void salvage(TTrack &track, AList< TMLink > &list) const
salvages links in a list. Used links will be removed from a list.
Definition TBuilder0.cxx:1511
TBuilder0::fit
virtual int fit(TTrackBase &) const
fits a track using a private fitter.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TBuilder0.h:128
TBuilder0::appendClusters
void appendClusters(TTrack &track, const AList< TMLink > &) const
appends TMLinks in a list.
Definition TBuilder0.cxx:1460
TCircle
A class to represent a circle in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TCircle.h:42
TLine0
A class to represent a track in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TLine0.h:30
TLine0::fit2sp
int fit2sp()
Definition TLine0.cxx:330
TLine0::chi2
double chi2(void) const
returns chi2.
Definition TLine0.cxx:111
TLine0::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/TLine0.h:142
TLine0::appendByszdistance
void appendByszdistance(AList< TMLink > &list, unsigned isl, float maxSigma)
Definition TLine0.cxx:577
TLine0::b
double b(void) const
returns coefficient b.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TLine0.h:135
TLine0::reducedChi2
double reducedChi2(void) const
returns reduced-chi2.
Definition TLine0.cxx:618
TLine0::removeChits
void removeChits()
remove extremly bad points.
Definition TLine0.cxx:524
TLine0::removeSLY
void removeSLY(AList< TMLink > &list)
Definition TLine0.cxx:565
TLine0::fit2p
int fit2p()
fits itself using isolated hits. Error was happened if return value is not zero.
Definition TLine0.cxx:410
TLine0::refine
void refine(AList< TMLink > &list, float maxSigma)
Definition TLine0.cxx:132
TLine0::a
double a(void) const
returns coefficient a.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TLine0.h:128
TMDCLayer::axialStereoLayerId
unsigned axialStereoLayerId(void) const
returns id of axial or stereo id.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCLayer.h:132
TMDCLayer::nWires
unsigned nWires(void) const
returns # of wires.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCLayer.h:116
TMDCWireHit::drift
float drift(unsigned) const
returns drift distance.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:218
TMDCWireHit::wire
const TMDCWire *const wire(void) const
returns a pointer to a TMDCWire.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:212
TMDCWireHit::xyPosition
const HepPoint3D & xyPosition(void) const
returns drift time
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:234
TMDCWire::consective
bool consective(const TMDCWire &) const
returns true if a given wire is consective in a layer.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:296
TMDCWire::layer
const TMDCLayer *const layer(void) const
returns a pointer to a layer.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:215
TMDCWire::localId
unsigned localId(void) const
returns local id in a wire layer.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:205
TMDCWire::layerId
unsigned layerId(void) const
returns layer id.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:207
TMDCWire::superLayerId
unsigned superLayerId(void) const
returns super layer id.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:209
TMDCWire::name
std::string name(void) const
returns name.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:308
TMDC::getTMDC
static TMDC * getTMDC(void)
Definition TMDC.cxx:84
TMDC::superLayer
const AList< TMDCLayer > *const superLayer(unsigned id) const
returns a pointer to a super-layer. 0 will be returned if 'id' is invalid.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDC.h:252
TMSelector::maxDistance
double maxDistance(void) const
returns max. distance required for stereo hits.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMSelector.h:205
TMSelector::nSuperLayers
unsigned nSuperLayers(void) const
returns min. # of super layers required.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMSelector.h:141
TMSelector::maxImpact
double maxImpact(void) const
returns max. impact(2D) required.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMSelector.h:166
TMSelector::nLinks
unsigned nLinks(void) const
returns min. # of hits(TMLinks) requried.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMSelector.h:129
TMSelector::maxSigma
double maxSigma(void) const
returns max. sigma for each TMLink.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMSelector.h:179
TMSelector::nLinksStereo
unsigned nLinksStereo(void) const
returns min. # of stereo hits(TMLinks) requried.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMSelector.h:192
TMSelector::minPt
double minPt(void) const
returns min. pt required.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMSelector.h:154
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::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::appendByApproach
void appendByApproach(AList< TMLink > &list, double maxSigma)
Definition TTrackBase.cxx:96
TTrack
A class to represent a track in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:128
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
EvtCyclic3::next
Index next(Index i)
Definition EvtCyclic3.cc:105
m2
double double * m2
Definition qcdloop1.h:83
t
int t()
Definition t.c:1