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