00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "TrkReco/TLine0.h"
00014 #include "TrkReco/TMDCUtil.h"
00015 #include "TrkReco/TMDCWire.h"
00016 #include "TrkReco/TMDCWireHit.h"
00017 #include "TrkReco/TMDCWireHitMC.h"
00018 #include "TrkReco/TTrackHEP.h"
00019
00020 const TLineFitter
00021 TLine0::_fitter = TLineFitter("TLine0 Default Line Fitter");
00022
00023 TLine0::TLine0()
00024 : TTrackBase(),
00025 _a(0.),
00026 _b(0.),
00027 _det(0.),
00028 _fittedUpdated(false),
00029 _chi2(0.),
00030 _reducedChi2(0.) {
00031
00032
00033 fitter(& TLine0::_fitter);
00034 }
00035
00036 TLine0::TLine0(const AList<TMLink> & a)
00037 : TTrackBase(a),
00038 _a(0.),
00039 _b(0.),
00040 _det(0.),
00041 _fittedUpdated(false),
00042 _chi2(0.),
00043 _reducedChi2(0.) {
00044
00045
00046 fitter(& TLine0::_fitter);
00047 }
00048
00049 TLine0::~TLine0() {
00050 }
00051
00052 void
00053 TLine0::dump(const std::string & msg, const std::string & pre) const {
00054 bool def = false;
00055 if (msg == "") def = true;
00056
00057 if (def || msg.find("line") != std::string::npos || msg.find("detail") != std::string::npos) {
00058 std::cout << pre;
00059 std::cout << "#links=" << _links.length();
00060 std::cout << ",a=" << _a;
00061 std::cout << ",b=" << _b;
00062 std::cout << ",det=" << _det;
00063 std::cout << std::endl;
00064 }
00065 if (! def) TTrackBase::dump(msg, pre);
00066 }
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 double
00113 TLine0::chi2(void) const {
00114 #ifdef TRKRECO_DEBUG_DETAIL
00115 if (! _fitted) std::cout << "TLine0::chi2 !!! fit not performed" << std::endl;
00116 #endif
00117
00118 if (_fittedUpdated) return _chi2;
00119 _chi2 = 0.;
00120 unsigned n = _links.length();
00121 for (unsigned i = 0; i < n; i++) {
00122 TMLink & l = * _links[i];
00123
00124 double x = l.position().x();
00125 double y = l.position().y();
00126 double c = y - _a * x - _b;
00127 _chi2 += c * c;
00128 }
00129 _fittedUpdated = true;
00130 return _chi2;
00131 }
00132
00133 void
00134 TLine0::refine(AList<TMLink> & list, float maxSigma) {
00135 AList<TMLink> bad;
00136 unsigned n = _links.length();
00137 for (unsigned i = 0; i < n; i++) {
00138 TMLink & l = * _links[i];
00139 double dist = distance(l);
00140 if (dist > maxSigma) bad.append(l);
00141 }
00142
00143 #ifdef TRKRECO_DEBUG_DETAIL
00144 std::cout << " TLine0::refine ... rejected hits:max distance=" << maxSigma;
00145 std::cout << std::endl;
00146 bad.sort(SortByWireId);
00147 for (unsigned i = 0; i < bad.length(); i++) {
00148 TMLink & l = * _links[i];
00149 std::cout << " ";
00150 std::cout << l.wire()->layerId() << "-";
00151 std::cout << l.wire()->localId();
00152 std::cout << "(";
00153 if (l.hit()->mc()) {
00154 if (l.hit()->mc()->hep()) std::cout << l.hit()->mc()->hep()->id();
00155 else std::cout << "0";
00156 }
00157 std::cout << "),";
00158 std::cout << l.position() << "," << distance(l);
00159 if (distance(l) > maxSigma) std::cout << " X";
00160 std::cout << std::endl;
00161 }
00162 #endif
00163
00164 if (bad.length()) {
00165 _links.remove(bad);
00166 list.append(bad);
00167 _fitted = false;
00168 _fittedUpdated = false;
00169 }
00170 }
00171
00172 int
00173 TLine0::fit2() {
00174
00175
00176 unsigned n = _links.length();
00177 int mask[100] = {0};
00178 int nsl[11] = {64,80,96,128,144,160,176,208,240,256,288};
00179 int npos = 0, nneg = 0;
00180 for (unsigned i = 0; i < n -1 ; i++) {
00181 TMLink & l = * _links[i];
00182 for (unsigned j = i+1; j < n ; j++) {
00183 TMLink & s = * _links[j];
00184 if(l.hit()->wire()->layerId() == s.hit()->wire()->layerId()){
00185
00186 if(i > 0 && (mask[i-1] == 1 && mask[j] == 1) ){
00187 TMLink & t = * _links[i-1];
00188 if(l.hit()->wire()->layerId() == t.hit()->wire()->layerId()){
00189 mask[i] = 1;
00190 }
00191 }
00192 int ilast = nsl[l.hit()->wire()->superLayerId()]-1;
00193 int ilocal = l.hit()->wire()->localId();
00194 int jlocal = s.hit()->wire()->localId();
00195 if(ilocal > 0 && ilocal < ilast){
00196 if(abs(jlocal-ilocal) > 1 ) {
00197 mask[i] = 1;
00198 mask[j] = 1;
00199 }
00200 }else if(ilocal == 0){
00201 if(jlocal > 1 && jlocal < ilast){
00202 mask[i] = 1;
00203 mask[j] = 1;
00204 }
00205 }else if(ilocal == ilast){
00206 if(jlocal > 0 && jlocal < ilast-1){
00207 mask[i] = 1;
00208 mask[j] = 1;
00209 }
00210 }
00211 }
00212 }
00213
00214 if(mask[i] == 0){
00215 if(l.position().y() >= 0) npos += 1;
00216 if(l.position().y() < 0) nneg += 1;
00217 }
00218 }
00219
00220
00221 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
00222 int nused = 0;
00223 int lyid[2];
00224 for (unsigned i = 0; i < n; i++) {
00225
00226 if(mask[i] == 1) continue;
00227
00228 TMLink & l = * _links[i];
00229
00230 double x = l.position().x();
00231 double y = l.position().y();
00232 if(abs(npos-nneg) > 3){
00233 if(npos > nneg && y < 0 ) continue;
00234 if(npos < nneg && y > 0 ) continue;
00235 }
00236 sumX += x;
00237 sumY += y;
00238 sumX2 += x * x;
00239 sumXY += x * y;
00240 sumY2 += y * y;
00241 if(nused < 2){
00242 lyid[nused] = l.hit()->wire()->layerId();
00243 }
00244 nused += 1;
00245 }
00246
00247 if(nused < 2 || (nused == 2 && lyid[0]==lyid[1])) {
00248 return -2;
00249 }
00250 double sum = double(nused);
00251 _det = sum * sumX2 - sumX * sumX;
00252 if(_det == 0.) {
00253 return -1;
00254 }
00255 _a = (sumXY * sum - sumX * sumY) / _det;
00256 _b = (sumX2 * sumY - sumX * sumXY) / _det;
00257
00258 _fitted = true;
00259 return 0;
00260 }
00261
00262 int
00263 TLine0::fit2s() {
00264
00265
00266 unsigned n = _links.length();
00267 int mask[100] = {0};
00268 int npos = 0, nneg = 0;
00269 for (unsigned i = 0; i < n -1 ; i++) {
00270 TMLink & l = * _links[i];
00271 for (unsigned j = i+1; j < n ; j++) {
00272 TMLink & s = * _links[j];
00273 if(l.hit()->wire()->layerId() == s.hit()->wire()->layerId()){
00274 mask[i] = 1;
00275 mask[j] = 1;
00276 }
00277 }
00278
00279 if(mask[i] == 0){
00280 if(l.position().y() >= 0) npos += 1;
00281 if(l.position().y() < 0) nneg += 1;
00282 }
00283 }
00284
00285
00286 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
00287 int nused = 0;
00288 int lyid[2];
00289 for (unsigned i = 0; i < n; i++) {
00290
00291 if(mask[i] == 1) continue;
00292
00293 TMLink & l = * _links[i];
00294
00295 double x = l.position().x();
00296 double y = l.position().y();
00297 if(npos > nneg && y < 0 ) continue;
00298 if(npos < nneg && y > 0 ) continue;
00299
00300 sumX += x;
00301 sumY += y;
00302 sumX2 += x * x;
00303 sumXY += x * y;
00304 sumY2 += y * y;
00305 nused += 1;
00306 }
00307
00308 if(nused < 4) {
00309 return -2;
00310 }
00311 double sum = double(nused);
00312 _det = sum * sumX2 - sumX * sumX;
00313 if(_det == 0.) {
00314 return -1;
00315 }
00316 _a = (sumXY * sum - sumX * sumY) / _det;
00317 _b = (sumX2 * sumY - sumX * sumXY) / _det;
00318
00319 _fitted = true;
00320 return 0;
00321 }
00322 int
00323 TLine0::fit2sp() {
00324
00325
00326 unsigned n = _links.length();
00327 int mask[100] = {0};
00328 double phi_ave = 0.;
00329 int nphi = 0;
00330 double Crad = 180./3.141592;
00331 for (unsigned i = 0; i < n -1 ; i++) {
00332 TMLink & l = * _links[i];
00333 for (unsigned j = i+1; j < n ; j++) {
00334 TMLink & s = * _links[j];
00335 if(l.hit()->wire()->layerId() == s.hit()->wire()->layerId()){
00336 mask[i] = 1;
00337 mask[j] = 1;
00338 }
00339 }
00340
00341 if(mask[i] != 1){
00342 double phi = Crad*atan2(l.position().y(), l.position().x());
00343 phi_ave += phi;
00344 nphi += 1;
00345 }
00346 }
00347
00348
00349 if(mask[n-1] != 1){
00350 TMLink & l = * _links[n-1];
00351 double phi = Crad*atan2(l.position().y(), l.position().x());
00352 phi_ave += phi;
00353 nphi += 1;
00354 }
00355 double phi_max = 0.;
00356 double phi_min = 0.;
00357 if(nphi> 0){
00358 phi_max = phi_ave/n + 40;
00359 phi_min = phi_ave/n - 40;
00360 }
00361
00362
00363 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
00364 int nused = 0;
00365 int lyid[2];
00366 for (unsigned i = 0; i < n; i++) {
00367
00368 if(mask[i] == 1) continue;
00369
00370 TMLink & l = * _links[i];
00371
00372 double x = l.position().x();
00373 double y = l.position().y();
00374 double phi = Crad*atan2(l.position().y(), l.position().x());
00375 if(phi > phi_max && phi<phi_min ) continue;
00376
00377 sumX += x;
00378 sumY += y;
00379 sumX2 += x * x;
00380 sumXY += x * y;
00381 sumY2 += y * y;
00382 nused += 1;
00383 }
00384
00385 if(nused < 4) {
00386 return -2;
00387 }
00388 double sum = double(nused);
00389 _det = sum * sumX2 - sumX * sumX;
00390 if(_det == 0.) {
00391 return -1;
00392 }
00393 _a = (sumXY * sum - sumX * sumY) / _det;
00394 _b = (sumX2 * sumY - sumX * sumXY) / _det;
00395
00396 _fitted = true;
00397 return 0;
00398 }
00399
00400 int
00401 TLine0::fit2p() {
00402
00403
00404 unsigned n = _links.length();
00405 int mask[100] = {0};
00406 int nsl[11] = {64,80,96,128,144,160,176,208,240,256,288};
00407 double phi_ave = 0.;
00408 int nphi = 0;
00409 double Crad = 180./3.141592;
00410 for (unsigned i = 0; i < n -1 ; i++) {
00411 TMLink & l = * _links[i];
00412 for (unsigned j = i+1; j < n ; j++) {
00413 TMLink & s = * _links[j];
00414 if(l.hit()->wire()->layerId() == s.hit()->wire()->layerId()){
00415
00416 if(i > 0 && (mask[i-1] == 1 && mask[j] == 1) ){
00417 TMLink & t = * _links[i-1];
00418 if(l.hit()->wire()->layerId() == t.hit()->wire()->layerId()){
00419 mask[i] = 1;
00420 }
00421 }
00422 int ilast = nsl[l.hit()->wire()->superLayerId()]-1;
00423 int ilocal = l.hit()->wire()->localId();
00424 int jlocal = s.hit()->wire()->localId();
00425 if(ilocal > 0 && ilocal < ilast){
00426 if(abs(jlocal-ilocal) > 1 ) {
00427 mask[i] = 1;
00428 mask[j] = 1;
00429 }
00430 }else if(ilocal == 0){
00431 if(jlocal > 1 && jlocal < ilast){
00432 mask[i] = 1;
00433 mask[j] = 1;
00434 }
00435 }else if(ilocal == ilast){
00436 if(jlocal > 0 && jlocal < ilast-1){
00437 mask[i] = 1;
00438 mask[j] = 1;
00439 }
00440 }
00441 }
00442 }
00443
00444
00445 if(mask[i] != 1){
00446 double phi = Crad*atan2(l.position().y(), l.position().x());
00447 phi_ave += phi;
00448 nphi += 1;
00449 }
00450 }
00451
00452
00453 if(mask[n-1] != 1){
00454 TMLink & l = * _links[n-1];
00455 double phi = Crad*atan2(l.position().y(), l.position().x());
00456 phi_ave += phi;
00457 nphi += 1;
00458 }
00459 double phi_max = 0.;
00460 double phi_min = 0.;
00461 if(nphi> 0){
00462 phi_max = phi_ave/n + 40;
00463 phi_min = phi_ave/n - 40;
00464 }
00465
00466
00467 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
00468 int nused = 0;
00469 int lyid[2];
00470 for (unsigned i = 0; i < n; i++) {
00471
00472 if(mask[i] == 1) continue;
00473
00474 TMLink & l = * _links[i];
00475
00476 double x = l.position().x();
00477 double y = l.position().y();
00478 double phi = Crad*atan2(l.position().y(), l.position().x());
00479 if(phi > phi_max && phi<phi_min ) continue;
00480
00481 sumX += x;
00482 sumY += y;
00483 sumX2 += x * x;
00484 sumXY += x * y;
00485 sumY2 += y * y;
00486 if(nused < 2){
00487 lyid[nused] = l.hit()->wire()->layerId();
00488 }
00489 nused += 1;
00490 }
00491
00492 if(nused < 2 || (nused == 2 && lyid[0]==lyid[1])) {
00493 return -2;
00494 }
00495 double sum = double(nused);
00496 _det = sum * sumX2 - sumX * sumX;
00497 if(_det == 0.) {
00498 return -1;
00499 }
00500 _a = (sumXY * sum - sumX * sumY) / _det;
00501 _b = (sumX2 * sumY - sumX * sumXY) / _det;
00502
00503 _fitted = true;
00504 return 0;
00505 }
00506
00507 void
00508 TLine0::removeChits() {
00509
00510 unsigned n = _links.length();
00511 int nlyr[50] = {0};
00512 int nneg = 0, npos = 0;
00513 for (unsigned i = 0; i < n -1 ; i++) {
00514 TMLink & l = * _links[i];
00515 nlyr[l.hit()->wire()->layerId()] += 1;
00516 if(l.position().y() < 0.){
00517 nneg += 1;
00518 }else {
00519 npos += 1;
00520 }
00521 }
00522
00523
00524 AList<TMLink> bad;
00525 for (unsigned i = 0; i < n; i++) {
00526
00527 TMLink & l = * _links[i];
00528
00529
00530 if (nlyr[l.hit()->wire()->layerId()] > 3) {
00531 bad.append(l);
00532 continue;
00533 }
00534
00535 if(abs(nneg - npos) > 3) {
00536 if(npos > nneg && l.position().y() < 0 ) bad.append(l);
00537 if(npos < nneg && l.position().y() > 0 ) bad.append(l);
00538 }
00539 }
00540
00541 if (bad.length() > 0 && bad.length() < n) {
00542 _links.remove(bad);
00543 }
00544
00545
00546 _fitted = false;
00547 _fittedUpdated = false;
00548 }
00549
00550 void
00551 TLine0::removeSLY(AList<TMLink> & list) {
00552 _links.remove(list);
00553 _fitted = false;
00554 _fittedUpdated = false;
00555 }
00556
00557 void
00558 TLine0::appendSLY(AList<TMLink> & list) {
00559 _links.append(list);
00560 _fitted = false;
00561 _fittedUpdated = false;
00562 }
00563
00564 void
00565 TLine0::appendByszdistance(AList<TMLink> & list, unsigned isl, float maxSigma) {
00566
00567
00568 unsigned nb = _links.length();
00569
00570
00571 unsigned n = list.length();
00572 for (unsigned i = 0; i < n; i++) {
00573 TMLink & l = * list[i];
00574 if(l.hit()->wire()->superLayerId() == isl){
00575 double dist = distance(l);
00576 if (dist < maxSigma) {
00577 _links.append(l);
00578 }
00579 }
00580 }
00581
00582 unsigned na = _links.length();
00583 if(nb != na){
00584 AList<TMLink> bad;
00585
00586 for(unsigned i = 0; i<na ;i++){
00587 TMLink & l = * _links[i];
00588 if(i < na - 1) {
00589 TMLink & lnext = * _links[i+1];
00590 if(l.hit()->wire()->layerId() == lnext.hit()->wire()->layerId() ){
00591 if(l.hit()->wire()->localId() == lnext.hit()->wire()->localId()){
00592 bad.append(l);
00593 }
00594 }
00595 }
00596 }
00597 if(bad.length() > 0) _links.remove(bad);
00598 _fitted = false;
00599 _fittedUpdated = false;
00600 }
00601 }
00602
00603 double
00604 TLine0::reducedChi2(void) const {
00605 #ifdef TRKRECO_DEBUG_DETAIL
00606 if (! _fitted) std::cout << "TLine0::reducedChi2 !!! fit not performed" << std::endl;
00607 #endif
00608
00609 if (_fittedUpdated) return _reducedChi2;
00610 double chi2 = 0.;
00611 double scale = 20.;
00612 unsigned n = _links.length();
00613 for (unsigned i = 0; i < n; i++) {
00614 TMLink & l = * _links[i];
00615
00616 double x = l.position().x();
00617 double y = l.position().y();
00618 double c = y - _a * x - _b;
00619 double err = scale * l.hit()->dDrift();
00620 chi2 += c * c / err / err;
00621 }
00622
00623 _reducedChi2 = chi2/(n-2);
00624 _fittedUpdated = true;
00625 return _reducedChi2;
00626 }