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