BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
TLine0.cxx
Go to the documentation of this file.
1//-----------------------------------------------------------------------------
2// $Id: TLine0.cxx,v 1.6 2010/03/31 09:58:59 liucy Exp $
3//-----------------------------------------------------------------------------
4// Filename : TLine0.cc
5// Section : Tracking
6// Owner : Yoshi Iwasaki
7// Email : yoshihito.iwasaki@kek.jp
8//-----------------------------------------------------------------------------
9// Description : A class to represent a line in tracking.
10// See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
11//-----------------------------------------------------------------------------
12
13#include "TrkReco/TLine0.h"
14#include "TrkReco/TMDCUtil.h"
15#include "TrkReco/TMDCWire.h"
16#include "TrkReco/TMDCWireHit.h"
17#include "TrkReco/TMDCWireHitMC.h"
18#include "TrkReco/TTrackHEP.h"
19
20const TLineFitter TLine0::_fitter = TLineFitter( "TLine0 Default Line Fitter" );
21
23 : TTrackBase()
24 , _a( 0. )
25 , _b( 0. )
26 , _det( 0. )
27 , _fittedUpdated( false )
28 , _chi2( 0. )
29 , _reducedChi2( 0. ) {
30
31 //...Set a defualt fitter...
32 fitter( &TLine0::_fitter );
33}
34
36 : TTrackBase( a )
37 , _a( 0. )
38 , _b( 0. )
39 , _det( 0. )
40 , _fittedUpdated( false )
41 , _chi2( 0. )
42 , _reducedChi2( 0. ) {
43
44 //...Set a defualt fitter...
45 fitter( &TLine0::_fitter );
46}
47
49
50void TLine0::dump( const std::string& msg, const std::string& pre ) const {
51 bool def = false;
52 if ( msg == "" ) def = true;
53
54 if ( def || msg.find( "line" ) != std::string::npos ||
55 msg.find( "detail" ) != std::string::npos )
56 {
57 std::cout << pre;
58 std::cout << "#links=" << _links.length();
59 std::cout << ",a=" << _a;
60 std::cout << ",b=" << _b;
61 std::cout << ",det=" << _det;
62 std::cout << std::endl;
63 }
64 if ( !def ) TTrackBase::dump( msg, pre );
65}
66
67// int
68// TLine0::fitx(void) {
69// if (_fitted) return 0;
70
71// unsigned n = _links.length();
72// double sum = double(n);
73// double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
74// for (unsigned i = 0; i < n; i++) {
75// TMLink & l = * _links[i];
76
77// double x = l.position().x();
78// double y = l.position().y();
79// sumX += x;
80// sumY += y;
81// sumX2 += x * x;
82// sumXY += x * y;
83// sumY2 += y * y;
84// }
85
86// _det = sum * sumX2 - sumX * sumX;
87// #ifdef TRKRECO_DEBUG_DETAIL
88// std::cout << " TLine0::fit ... det=" << _det << std::endl;
89// #endif
90// if(_det == 0. && n != 2){
91// return -1;
92// }else if(_det == 0. && n == 2){
93// double x0 = _links[0]->position().x();
94// double y0 = _links[0]->position().y();
95// double x1 = _links[1]->position().x();
96// double y1 = _links[1]->position().y();
97// if(x0 == x1)return -1;
98// _a = (y0-y1)/(x0-x1);
99// _b = -_a*x1 + y1;
100
101// _fitted = true;
102// return 0;
103// }
104// _a = (sumXY * sum - sumX * sumY) / _det;
105// _b = (sumX2 * sumY - sumX * sumXY) / _det;
106
107// _fitted = true;
108// return 0;
109// }
110
111double TLine0::chi2( void ) const {
112#ifdef TRKRECO_DEBUG_DETAIL
113 if ( !_fitted ) std::cout << "TLine0::chi2 !!! fit not performed" << std::endl;
114#endif
115
116 if ( _fittedUpdated ) return _chi2;
117 _chi2 = 0.;
118 unsigned n = _links.length();
119 for ( unsigned i = 0; i < n; i++ )
120 {
121 TMLink& l = *_links[i];
122
123 double x = l.position().x();
124 double y = l.position().y();
125 double c = y - _a * x - _b;
126 _chi2 += c * c;
127 }
128 _fittedUpdated = true;
129 return _chi2;
130}
131
132void TLine0::refine( AList<TMLink>& list, float maxSigma ) {
133 AList<TMLink> bad;
134 unsigned n = _links.length();
135 for ( unsigned i = 0; i < n; i++ )
136 {
137 TMLink& l = *_links[i];
138 double dist = distance( l );
139 if ( dist > maxSigma ) bad.append( l );
140 }
141
142#ifdef TRKRECO_DEBUG_DETAIL
143 std::cout << " TLine0::refine ... rejected hits:max distance=" << maxSigma;
144 std::cout << std::endl;
145 bad.sort( SortByWireId );
146 for ( unsigned i = 0; i < bad.length(); i++ )
147 {
148 TMLink& l = *_links[i];
149 std::cout << " ";
150 std::cout << l.wire()->layerId() << "-";
151 std::cout << l.wire()->localId();
152 std::cout << "(";
153 if ( l.hit()->mc() )
154 {
155 if ( l.hit()->mc()->hep() ) std::cout << l.hit()->mc()->hep()->id();
156 else std::cout << "0";
157 }
158 std::cout << "),";
159 std::cout << l.position() << "," << distance( l );
160 if ( distance( l ) > maxSigma ) std::cout << " X";
161 std::cout << std::endl;
162 }
163#endif
164
165 if ( bad.length() )
166 {
167 _links.remove( bad );
168 list.append( bad );
169 _fitted = false;
170 _fittedUpdated = false;
171 }
172}
173
175 // if (_fitted) return 0;
176
177 unsigned n = _links.length();
178 int mask[100] = { 0 };
179 int nsl[11] = { 64, 80, 96, 128, 144, 160, 176, 208, 240, 256, 288 };
180 int npos = 0, nneg = 0;
181 for ( unsigned i = 0; i < n - 1; i++ )
182 {
183 TMLink& l = *_links[i];
184 for ( unsigned j = i + 1; j < n; j++ )
185 {
186 TMLink& s = *_links[j];
187 if ( l.hit()->wire()->layerId() == s.hit()->wire()->layerId() )
188 {
189 //... Check 3 consective hits
190 if ( i > 0 && ( mask[i - 1] == 1 && mask[j] == 1 ) )
191 {
192 TMLink& t = *_links[i - 1];
193 if ( l.hit()->wire()->layerId() == t.hit()->wire()->layerId() ) { mask[i] = 1; }
194 }
195 int ilast = nsl[l.hit()->wire()->superLayerId()] - 1;
196 int ilocal = l.hit()->wire()->localId();
197 int jlocal = s.hit()->wire()->localId();
198 if ( ilocal > 0 && ilocal < ilast )
199 {
200 if ( abs( jlocal - ilocal ) > 1 )
201 {
202 mask[i] = 1;
203 mask[j] = 1;
204 }
205 }
206 else if ( ilocal == 0 )
207 {
208 if ( jlocal > 1 && jlocal < ilast )
209 {
210 mask[i] = 1;
211 mask[j] = 1;
212 }
213 }
214 else if ( ilocal == ilast )
215 {
216 if ( jlocal > 0 && jlocal < ilast - 1 )
217 {
218 mask[i] = 1;
219 mask[j] = 1;
220 }
221 }
222 }
223 }
224 //...
225 if ( mask[i] == 0 )
226 {
227 if ( l.position().y() >= 0 ) npos += 1;
228 if ( l.position().y() < 0 ) nneg += 1;
229 }
230 }
231
232 //....
233 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
234 int nused = 0;
235 int lyid[2];
236 for ( unsigned i = 0; i < n; i++ )
237 {
238
239 if ( mask[i] == 1 ) continue;
240
241 TMLink& l = *_links[i];
242
243 double x = l.position().x();
244 double y = l.position().y();
245 if ( abs( npos - nneg ) > 3 )
246 {
247 if ( npos > nneg && y < 0 ) continue;
248 if ( npos < nneg && y > 0 ) continue;
249 }
250 sumX += x;
251 sumY += y;
252 sumX2 += x * x;
253 sumXY += x * y;
254 sumY2 += y * y;
255 if ( nused < 2 ) { lyid[nused] = l.hit()->wire()->layerId(); }
256 nused += 1;
257 }
258
259 if ( nused < 2 || ( nused == 2 && lyid[0] == lyid[1] ) ) { return -2; }
260 double sum = double( nused );
261 _det = sum * sumX2 - sumX * sumX;
262 if ( _det == 0. ) { return -1; }
263 _a = ( sumXY * sum - sumX * sumY ) / _det;
264 _b = ( sumX2 * sumY - sumX * sumXY ) / _det;
265
266 _fitted = true;
267 return 0;
268}
269
271 // if (_fitted) return 0;
272
273 unsigned n = _links.length();
274 int mask[100] = { 0 };
275 int npos = 0, nneg = 0;
276 for ( unsigned i = 0; i < n - 1; i++ )
277 {
278 TMLink& l = *_links[i];
279 for ( unsigned j = i + 1; j < n; j++ )
280 {
281 TMLink& s = *_links[j];
282 if ( l.hit()->wire()->layerId() == s.hit()->wire()->layerId() )
283 {
284 mask[i] = 1;
285 mask[j] = 1;
286 }
287 }
288 //...
289 if ( mask[i] == 0 )
290 {
291 if ( l.position().y() >= 0 ) npos += 1;
292 if ( l.position().y() < 0 ) nneg += 1;
293 }
294 }
295
296 //....
297 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
298 int nused = 0;
299 int lyid[2];
300 for ( unsigned i = 0; i < n; i++ )
301 {
302
303 if ( mask[i] == 1 ) continue;
304
305 TMLink& l = *_links[i];
306
307 double x = l.position().x();
308 double y = l.position().y();
309 if ( npos > nneg && y < 0 ) continue;
310 if ( npos < nneg && y > 0 ) continue;
311
312 sumX += x;
313 sumY += y;
314 sumX2 += x * x;
315 sumXY += x * y;
316 sumY2 += y * y;
317 nused += 1;
318 }
319
320 if ( nused < 4 ) { return -2; }
321 double sum = double( nused );
322 _det = sum * sumX2 - sumX * sumX;
323 if ( _det == 0. ) { return -1; }
324 _a = ( sumXY * sum - sumX * sumY ) / _det;
325 _b = ( sumX2 * sumY - sumX * sumXY ) / _det;
326
327 _fitted = true;
328 return 0;
329}
331 // if (_fitted) return 0;
332
333 unsigned n = _links.length();
334 int mask[100] = { 0 };
335 double phi_ave = 0.;
336 int nphi = 0;
337 double Crad = 180. / 3.141592;
338 for ( unsigned i = 0; i < n - 1; i++ )
339 {
340 TMLink& l = *_links[i];
341 for ( unsigned j = i + 1; j < n; j++ )
342 {
343 TMLink& s = *_links[j];
344 if ( l.hit()->wire()->layerId() == s.hit()->wire()->layerId() )
345 {
346 mask[i] = 1;
347 mask[j] = 1;
348 }
349 }
350 //...
351 if ( mask[i] != 1 )
352 {
353 double phi = Crad * atan2( l.position().y(), l.position().x() );
354 phi_ave += phi;
355 nphi += 1;
356 }
357 }
358
359 //...
360 if ( mask[n - 1] != 1 )
361 {
362 TMLink& l = *_links[n - 1];
363 double phi = Crad * atan2( l.position().y(), l.position().x() );
364 phi_ave += phi;
365 nphi += 1;
366 }
367 double phi_max = 0.;
368 double phi_min = 0.;
369 if ( nphi > 0 )
370 {
371 phi_max = phi_ave / n + 40;
372 phi_min = phi_ave / n - 40;
373 }
374
375 //....
376 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
377 int nused = 0;
378 int lyid[2];
379 for ( unsigned i = 0; i < n; i++ )
380 {
381
382 if ( mask[i] == 1 ) continue;
383
384 TMLink& l = *_links[i];
385
386 double x = l.position().x();
387 double y = l.position().y();
388 double phi = Crad * atan2( l.position().y(), l.position().x() );
389 if ( phi > phi_max && phi < phi_min ) continue;
390
391 sumX += x;
392 sumY += y;
393 sumX2 += x * x;
394 sumXY += x * y;
395 sumY2 += y * y;
396 nused += 1;
397 }
398
399 if ( nused < 4 ) { return -2; }
400 double sum = double( nused );
401 _det = sum * sumX2 - sumX * sumX;
402 if ( _det == 0. ) { return -1; }
403 _a = ( sumXY * sum - sumX * sumY ) / _det;
404 _b = ( sumX2 * sumY - sumX * sumXY ) / _det;
405
406 _fitted = true;
407 return 0;
408}
409
411 // if (_fitted) return 0;
412
413 unsigned n = _links.length();
414 int mask[100] = { 0 };
415 int nsl[11] = { 64, 80, 96, 128, 144, 160, 176, 208, 240, 256, 288 };
416 double phi_ave = 0.;
417 int nphi = 0;
418 double Crad = 180. / 3.141592;
419 for ( unsigned i = 0; i < n - 1; i++ )
420 {
421 TMLink& l = *_links[i];
422 for ( unsigned j = i + 1; j < n; j++ )
423 {
424 TMLink& s = *_links[j];
425 if ( l.hit()->wire()->layerId() == s.hit()->wire()->layerId() )
426 {
427 //... Check 3 consective hits
428 if ( i > 0 && ( mask[i - 1] == 1 && mask[j] == 1 ) )
429 {
430 TMLink& t = *_links[i - 1];
431 if ( l.hit()->wire()->layerId() == t.hit()->wire()->layerId() ) { mask[i] = 1; }
432 }
433 int ilast = nsl[l.hit()->wire()->superLayerId()] - 1;
434 int ilocal = l.hit()->wire()->localId();
435 int jlocal = s.hit()->wire()->localId();
436 if ( ilocal > 0 && ilocal < ilast )
437 {
438 if ( abs( jlocal - ilocal ) > 1 )
439 {
440 mask[i] = 1;
441 mask[j] = 1;
442 }
443 }
444 else if ( ilocal == 0 )
445 {
446 if ( jlocal > 1 && jlocal < ilast )
447 {
448 mask[i] = 1;
449 mask[j] = 1;
450 }
451 }
452 else if ( ilocal == ilast )
453 {
454 if ( jlocal > 0 && jlocal < ilast - 1 )
455 {
456 mask[i] = 1;
457 mask[j] = 1;
458 }
459 }
460 }
461 }
462 //...
463 //...
464 if ( mask[i] != 1 )
465 {
466 double phi = Crad * atan2( l.position().y(), l.position().x() );
467 phi_ave += phi;
468 nphi += 1;
469 }
470 }
471
472 //...
473 if ( mask[n - 1] != 1 )
474 {
475 TMLink& l = *_links[n - 1];
476 double phi = Crad * atan2( l.position().y(), l.position().x() );
477 phi_ave += phi;
478 nphi += 1;
479 }
480 double phi_max = 0.;
481 double phi_min = 0.;
482 if ( nphi > 0 )
483 {
484 phi_max = phi_ave / n + 40;
485 phi_min = phi_ave / n - 40;
486 }
487
488 //....
489 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
490 int nused = 0;
491 int lyid[2];
492 for ( unsigned i = 0; i < n; i++ )
493 {
494
495 if ( mask[i] == 1 ) continue;
496
497 TMLink& l = *_links[i];
498
499 double x = l.position().x();
500 double y = l.position().y();
501 double phi = Crad * atan2( l.position().y(), l.position().x() );
502 if ( phi > phi_max && phi < phi_min ) continue;
503
504 sumX += x;
505 sumY += y;
506 sumX2 += x * x;
507 sumXY += x * y;
508 sumY2 += y * y;
509 if ( nused < 2 ) { lyid[nused] = l.hit()->wire()->layerId(); }
510 nused += 1;
511 }
512
513 if ( nused < 2 || ( nused == 2 && lyid[0] == lyid[1] ) ) { return -2; }
514 double sum = double( nused );
515 _det = sum * sumX2 - sumX * sumX;
516 if ( _det == 0. ) { return -1; }
517 _a = ( sumXY * sum - sumX * sumY ) / _det;
518 _b = ( sumX2 * sumY - sumX * sumXY ) / _det;
519
520 _fitted = true;
521 return 0;
522}
523
525
526 unsigned n = _links.length();
527 int nlyr[50] = { 0 };
528 int nneg = 0, npos = 0;
529 for ( unsigned i = 0; i < n - 1; i++ )
530 {
531 TMLink& l = *_links[i];
532 nlyr[l.hit()->wire()->layerId()] += 1;
533 if ( l.position().y() < 0. ) { nneg += 1; }
534 else { npos += 1; }
535 }
536
537 //...
538 AList<TMLink> bad;
539 for ( unsigned i = 0; i < n; i++ )
540 {
541
542 TMLink& l = *_links[i];
543
544 //...if # of hits in a wire layer, don't use...
545 if ( nlyr[l.hit()->wire()->layerId()] > 3 )
546 {
547 bad.append( l );
548 continue;
549 }
550 //...remove extremely bad poinits
551 if ( abs( nneg - npos ) > 3 )
552 {
553 if ( npos > nneg && l.position().y() < 0 ) bad.append( l );
554 if ( npos < nneg && l.position().y() > 0 ) bad.append( l );
555 }
556 }
557 //...
558 if ( bad.length() > 0 && bad.length() < n ) { _links.remove( bad ); }
559
560 //... For the next fit
561 _fitted = false;
562 _fittedUpdated = false;
563}
564
566 _links.remove( list );
567 _fitted = false;
568 _fittedUpdated = false;
569}
570
572 _links.append( list );
573 _fitted = false;
574 _fittedUpdated = false;
575}
576
577void TLine0::appendByszdistance( AList<TMLink>& list, unsigned isl, float maxSigma ) {
578
579 //... intialize
580 unsigned nb = _links.length();
581
582 //....Select good hit
583 unsigned n = list.length();
584 for ( unsigned i = 0; i < n; i++ )
585 {
586 TMLink& l = *list[i];
587 if ( l.hit()->wire()->superLayerId() == isl )
588 {
589 double dist = distance( l );
590 if ( dist < maxSigma ) { _links.append( l ); }
591 }
592 }
593
594 unsigned na = _links.length();
595 if ( nb != na )
596 {
597 AList<TMLink> bad;
598 //... remove duplicated hits
599 for ( unsigned i = 0; i < na; i++ )
600 {
601 TMLink& l = *_links[i];
602 if ( i < na - 1 )
603 {
604 TMLink& lnext = *_links[i + 1];
605 if ( l.hit()->wire()->layerId() == lnext.hit()->wire()->layerId() )
606 {
607 if ( l.hit()->wire()->localId() == lnext.hit()->wire()->localId() )
608 { bad.append( l ); }
609 }
610 }
611 }
612 if ( bad.length() > 0 ) _links.remove( bad );
613 _fitted = false;
614 _fittedUpdated = false;
615 }
616}
617
618double TLine0::reducedChi2( void ) const {
619#ifdef TRKRECO_DEBUG_DETAIL
620 if ( !_fitted ) std::cout << "TLine0::reducedChi2 !!! fit not performed" << std::endl;
621#endif
622
623 if ( _fittedUpdated ) return _reducedChi2;
624 double chi2 = 0.;
625 double scale = 20.;
626 unsigned n = _links.length();
627 for ( unsigned i = 0; i < n; i++ )
628 {
629 TMLink& l = *_links[i];
630
631 double x = l.position().x();
632 double y = l.position().y();
633 double c = y - _a * x - _b;
634 double err = scale * l.hit()->dDrift();
635 chi2 += c * c / err / err;
636 }
637
638 _reducedChi2 = chi2 / ( n - 2 );
639 _fittedUpdated = true;
640 return _reducedChi2;
641}
int nneg
const Int_t n
XmlRpcServer s
void appendSLY(AList< TMLink > &list)
Definition TLine0.cxx:571
int fit2sp()
Definition TLine0.cxx:330
double chi2(void) const
returns chi2.
Definition TLine0.cxx:111
double distance(const TMLink &) const
returns distance to a position of TMLink itself. (not to a wire)
void appendByszdistance(AList< TMLink > &list, unsigned isl, float maxSigma)
Definition TLine0.cxx:577
double reducedChi2(void) const
returns reduced-chi2.
Definition TLine0.cxx:618
int fit2()
fits itself. Error was happened if return value is not zero.
Definition TLine0.cxx:174
virtual ~TLine0()
Destructor.
Definition TLine0.cxx:48
void removeChits()
remove extremly bad points.
Definition TLine0.cxx:524
void removeSLY(AList< TMLink > &list)
Definition TLine0.cxx:565
int fit2p()
fits itself using isolated hits. Error was happened if return value is not zero.
Definition TLine0.cxx:410
void refine(AList< TMLink > &list, float maxSigma)
Definition TLine0.cxx:132
int fit2s()
Definition TLine0.cxx:270
TLine0()
Constructor.
Definition TLine0.cxx:22
double a(void) const
returns coefficient a.
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TLine0.cxx:50
A class to fit a TTrackBase object to a line.
const TTrackHEP *const hep(void) const
returns a pointer to a GEN_HEPEVT.
float dDrift(unsigned) const
returns drift distance error.
const TMDCWire *const wire(void) const
returns a pointer to a TMDCWire.
const TMDCWireHitMC *const mc(void) const
returns a pointer to TMDCWireHitMC.
unsigned localId(void) const
returns local id in a wire layer.
unsigned layerId(void) const
returns layer id.
unsigned superLayerId(void) const
returns super layer id.
virtual void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
const TMFitter *const fitter(void) const
returns a pointer to a default fitter.
TTrackBase()
Constructor.
unsigned id(void) const
returns an id started from 0.
int t()
Definition t.c:1