00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
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
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384 #include <stdexcept>
00385 #include <cmath>
00386
00387
00388 #include "StThreeVector.hh"
00389 #include "StThreeVectorF.hh"
00390 #include "StThreeVectorD.hh"
00391 #include "StPhysicalHelixD.hh"
00392 #include "THelixTrack.h"
00393
00394 #include "StHit.h"
00395
00396
00397 #include "StiKalmanTrack.h"
00398 #include "StiKalmanTrackFinder.h"
00399 #include "StiToolkit.h"
00400 #include "StiDetectorContainer.h"
00401 #include "StiHit.h"
00402 #include "StiKalmanTrackNode.h"
00403 #include "StiKalmanTrack.h"
00404 #include "StiDetector.h"
00405 #include "StiDetectorGroups.h"
00406 #include "StiDetectorBuilder.h"
00407 #include "StiPlacement.h"
00408 #include "StiMaterial.h"
00409 #include "StDetectorDbMaker/StiHitErrorCalculator.h"
00410 #include "StPhysicalHelixD.hh"
00411 #include "StHelix.hh"
00412 #include "StDetectorDbMaker/StiKalmanTrackFitterParameters.h"
00413 #include "StDetectorDbMaker/StiKalmanTrackFinderParameters.h"
00414 #include "StiHitContainer.h"
00415 #include "StiTrackNodeHelper.h"
00416 #include "StiUtilities/StiDebug.h"
00417 #if ROOT_VERSION_CODE < 331013
00418 #include "TCL.h"
00419 #else
00420 #include "TCernLib.h"
00421 #endif
00422 #include "StMessMgr.h"
00423 ostream& operator<<(ostream&, const StiHit&);
00424
00425 Factory<StiKalmanTrackNode>* StiKalmanTrack::trackNodeFactory = 0;
00426 int StiKalmanTrack::mgMaxRefiter = 100;
00427 int StiKalmanTrack::_debug = 0;
00428 int debugCount=0;
00429
00430
00431 static StiTrackNodeHelper sTNH;
00432 static double diff(const StiNodePars &p1,const StiNodeErrs &e1
00433 ,const StiNodePars &p2,const StiNodeErrs &e2,int &igor);
00434
00435
00436
00437
00448
00449 void StiKalmanTrack::reset()
00450 {
00451 static int mIdCount = 0;
00452 if ((++mIdCount) >= 1<<16) mIdCount = 1;
00453 mIdDb = mIdCount;
00454 firstNode = 0;
00455 lastNode = 0;
00456 mSeedHitCount = 0;
00457 mVertex = 0;
00458 m = -1.;
00459 mFlag = 0;
00460 _dca = 0;
00461 _vChi2=-2;
00462 StiDebug::Break(mIdDb);
00463 }
00464
00465
00466
00467
00472
00473 void StiKalmanTrack::setKalmanTrackNodeFactory(Factory<StiKalmanTrackNode>* val)
00474 {
00475 trackNodeFactory = val;
00476 }
00477
00478
00479
00480
00512
00513 int StiKalmanTrack::initialize(const std::vector<StiHit*> &hits)
00514 {
00515 #ifdef DO_TPCCATRACKER
00516 initialize0(hits);
00517
00518 int ierr = approx(0);
00519 if (!ierr) return 0;
00520 BFactory::Free(this);
00521 return 1;
00522 }
00523
00524
00525 int StiKalmanTrack::initialize0(const std::vector<StiHit*> &hits, StiNodePars *firstPars, StiNodePars *lastPars, StiNodeErrs *firstErrs, StiNodeErrs *lastErrs)
00526 {
00527 #endif
00528
00529 reset();
00530
00531 const StiDetector* detector=0;
00532 UInt_t nhits = hits.size();
00533 setSeedHitCount(nhits);
00534 #ifdef DO_TPCCATRACKER
00535 StiDetectorContainer *detectorContainer = StiToolkit::instance()->getDetectorContainer();
00536 const StiDetector* detectorOld = 0;
00537 StiHit *hit_Old = 0;
00538 for (UInt_t ihit = 0; ihit < nhits; ihit++) {
00539 StiHit *hit = hits[ihit];
00540 detector = hit->detector();
00541 assert(detector);
00542
00543 if (hit_Old && detector->getGroupId() == kTpcId) {
00544 Double_t R_hit = detector->getPlacement()->getLayerRadius();
00545 Double_t angle_hit = detector->getPlacement()->getNormalRefAngle();
00546 detectorOld = hit_Old->detector();
00547 Double_t R_hit_OLD = detectorOld->getPlacement()->getLayerRadius();
00548 if (_debug && detectorOld == detector) {
00549 cout << "The same detector for hit " << ihit << endl;
00550 cout << "hit \t" << *hit << endl;
00551 if (hit_Old)
00552 cout << "hitOld\t" << *hit_Old << endl;
00553 }
00554 Double_t angle_hit_OLD = detectorOld->getPlacement()->getNormalRefAngle();
00555 if (TMath::Abs(angle_hit - angle_hit_OLD) < TMath::DegToRad()*5) {
00556 while ((R_hit < R_hit_OLD)) {
00557 detectorContainer->setToDetector( detectorOld );
00558 if ( detectorContainer->moveIn()) {
00559 StiDetector* d = detectorContainer->getCurrentDetector();
00560 if (d == detector) break;
00561 detectorOld = d;
00562 R_hit_OLD = detectorOld->getPlacement()->getLayerRadius();
00563 if (detectorOld->isActive()) {
00564 StiKalmanTrackNode * nI = trackNodeFactory->getInstance();
00565 nI->initialize(d);
00566 add(nI,kOutsideIn);
00567 }
00568 }
00569 }
00570 }
00571 }
00572 StiKalmanTrackNode * n = trackNodeFactory->getInstance();
00573 n->initialize(hit);
00574 add(n,kOutsideIn);
00575 detectorOld = (StiDetector*) detector;
00576 hit_Old = hit;
00577 }
00578 #else
00579 for (UInt_t ihit=0;ihit<nhits;ihit++)
00580 {
00581 StiHit *hit = hits[ihit];
00582 detector = hit->detector();
00583 assert(detector);
00584 StiKalmanTrackNode * n = trackNodeFactory->getInstance();
00585 n->initialize(hit);
00586 add(n,kOutsideIn);
00587 }
00588 #endif
00589 #ifdef DO_TPCCATRACKER
00590 if (firstPars){
00591 firstNode->fitPars() = *firstPars;
00592 }
00593 if (firstErrs){
00594 firstNode->fitErrs() = *firstErrs;
00595
00596 }
00597 if (lastPars){
00598 lastNode ->fitPars() = *lastPars;
00599 }
00600 if (lastErrs){
00601 lastNode->fitErrs() = *lastErrs;
00602
00603 }
00604 return 0;
00605 #else
00606 int ierr = approx(0);
00607 if (!ierr) return 0;
00608 BFactory::Free(this);
00609 return 1;
00610 #endif
00611 }
00612
00613
00614 StThreeVector<double> StiKalmanTrack::getMomentumAtOrigin() const
00615 {
00616 double px,py,pz;
00617 px=py=pz=0;
00618
00619 StiKalmanTrackNode * inner = getInnerMostNode();
00620
00621 if (inner==0)throw logic_error("StiKalmanTrack::getMomentumAtOrigin() - ERROR - No node");
00622 inner->propagate(0.,0,kOutsideIn);
00623 double p[3];
00624 inner->getMomentum(p,0);
00625 StThreeVector<double> p3(p[0],p[1],p[2]);
00626 p3.rotateZ(inner->getAlpha());
00627 return p3;
00628 }
00629
00630
00637
00638 int StiKalmanTrack::getCharge() const
00639 {
00640 StiKalmanTrackNode *node = getInnerMostNode();
00641 if (!node) return 0;
00642 return node->getCharge();
00643 }
00644
00645
00651 double StiKalmanTrack::getChi2() const
00652 {
00653 double fitHits = 0;
00654 double trackChi2 = 0;
00655 double nodeChi2 = 0;
00656 double maxChi2 = StiKalmanTrackFitterParameters::instance()->getMaxChi2();
00657 double theChi2 = 1.e+60;
00658 if (!firstNode) return theChi2;
00659 theChi2 = 0;
00660 StiKTNIterator it;
00661 for (it=begin();it!=end();it++) {
00662 StiKalmanTrackNode *node = &(*it);
00663 if (!node->isValid()) continue;
00664 if (!node->getHit() ) continue;
00665 nodeChi2 = node->getChi2();
00666 if (nodeChi2>maxChi2) continue;
00667 trackChi2 += nodeChi2;
00668 ++fitHits;
00669 }
00670 return (fitHits>5)?trackChi2/(2.*fitHits-5.):1e30;
00671 }
00672
00673
00683
00684 int StiKalmanTrack::getPointCount(int detectorId) const
00685 {
00686 const StiDetector *detector=0;
00687 int nPts = 0;
00688 StiKTNIterator it;
00689 for (it=begin();it!=end();it++) {
00690 StiKalmanTrackNode *node = &(*it);
00691 if (!node->isValid()) continue;
00692 if (!node->getHit()) continue;
00693 detector = node->getDetector();
00694 if (!detector) continue;
00695 if (detectorId && detector->getGroupId() != detectorId) continue;
00696 nPts++;
00697 }
00698 return nPts;
00699 }
00700
00701
00712
00713 int StiKalmanTrack::getMaxPointCount(int detectorId) const
00714 {
00715 int nPts = 0;
00716 StiKTNIterator it;
00717
00718 for (it=begin();it!=end();it++){
00719 const StiKalmanTrackNode *node = &(*it);
00720 if (!node->isValid()) continue;
00721 const StiDetector *detector = node->getDetector();
00722 if (!detector) continue;
00723 StiHit* h = node->getHit();
00724 if (!h && !detector->isActive(node->getY(),node->getZ())) continue;
00725 if (detectorId && detector->getGroupId() != detectorId) continue;
00726 nPts++;
00727 }
00728 return nPts;
00729 }
00730
00731
00732
00742 int StiKalmanTrack::getGapCount() const
00743 {
00744 int gaps = 0;
00745 if (firstNode)
00746 {
00747 StiKTNIterator it;
00748 bool inGap = false;
00749 for (it=begin();it!=end();it++)
00750 {
00751 const StiDetector * detector = (*it).getDetector();
00752 if (detector && detector->isActive())
00753 {
00754 if ((*it).getHit())
00755 {
00756 if (inGap)
00757 inGap = false;
00758 }
00759 else
00760 {
00761 if (!inGap)
00762 {
00763 inGap = true;
00764 gaps++;
00765 }
00766 }
00767 }
00768 }
00769 }
00770 return gaps;
00771 }
00772
00773
00783
00784
00785 int StiKalmanTrack::getFitPointCount(int detectorId) const
00786 {
00787 int fitPointCount = 0;
00788 StiKTNIterator it;
00789 for (it=begin();it!=end();it++) {
00790 StiKalmanTrackNode* node = &(*it);
00791 if(!node->isValid()) continue;
00792 StiHit* hit = node->getHit();
00793 if (!hit) continue;
00794 if (!node->isFitted()) continue;
00795 const StiDetector *det = hit->detector();
00796 if (!det) continue;
00797 if (detectorId && detectorId!=det->getGroupId())continue;
00798 fitPointCount++;
00799 }
00800 return fitPointCount;
00801 }
00802
00803 void StiKalmanTrack::getAllPointCount(int count[1][3],int maxDetId) const
00804 {
00805
00806
00807
00808
00809
00810
00811 enum {kPP=0,kMP=1,kFP=2};
00812
00813 memset(count[0],0,(maxDetId+1)*3*sizeof(int));
00814 StiKTNIterator it;
00815
00816 for (it=begin();it!=end();it++){
00817 const StiKalmanTrackNode *node = &(*it);
00818 if (!node->isValid()) continue;
00819 const StiDetector *detector = node->getDetector();
00820 if (!detector) continue;
00821 int detId = detector->getGroupId();
00822 StiHit* h = node->getHit();
00823
00824
00825 if (h || detector->isActive(node->getY(),node->getZ())) {
00826 count[0][kPP]++; count[detId][kPP]++;
00827 }
00828
00829 if (!h ) continue;
00830
00831 count[0][kMP]++; count[detId][kMP]++;
00832 if (!node->isFitted()) continue;
00833 count[0][kFP]++; count[detId][kFP]++;
00834 }
00835 }
00836
00837
00847 double StiKalmanTrack::getTrackLength() const
00848 {
00849 double x[2][4];
00850 for (int i = 0; i < 4; i++){ x[0][i] = 0; }
00851 double len=0;
00852 int iready=0;
00853 StiKalmanTrackNode *node;
00854 StiKTNIterator it = begin();
00855 for (;(node=it());it++){
00856 if (!node->isValid()) continue;
00857 if (!node->getHit()) continue;
00858 if ( node->getChi2()>10000.)continue;
00859 x[1][0]=node->x_g();
00860 x[1][1]=node->y_g();
00861 x[1][2]=node->z_g();
00862 x[1][3]=node->getCurvature();
00863 if (iready) {
00864 double dlen = sqrt(pow(x[1][0]-x[0][0],2) + pow(x[1][1]-x[0][1],2));
00865 double curv = fabs(0.5*(x[0][3]+x[1][3]));
00866 double dsin = (0.5*dlen*curv);
00867 if (dsin>0.9) {
00868 LOG_DEBUG <<
00869 Form("StiKalmanTrack::getTrackLength ***ERROR*** dsin %g >.9",dsin)
00870 << endm;
00871 dsin = 0.9;
00872 }
00873 dlen = (dsin<0.1)? dlen*(1.+dsin*dsin/6) : 2*asin(dsin)/curv;
00874 len +=sqrt(dlen*dlen + pow(x[1][2]-x[0][2],2));
00875 }
00876 memcpy(x[0],x[1],4*sizeof(double)); iready=2005;
00877 }
00878 return len;
00879 }
00880
00881
00889
00890
00891
00892
00893
00897
00898 double StiKalmanTrack::getTrackRadLength() const
00899 {
00900 double x1, x2, x3;
00901 double totalR=0.;
00902
00903
00904 StiKTNIterator tNode = begin();
00905
00906
00907
00908 StiKalmanTrackNode *thisNode = &(*tNode);
00909
00910 x1=thisNode->pathlength()/2.;
00911 x3=0.;
00912
00913
00914
00915 while ((tNode++)!=end() && (*tNode).getDetector())
00916 {
00917
00918 StiKalmanTrackNode *nextNode = &(*(tNode));
00919
00920
00921 x2=thisNode->pathLToNode(nextNode);
00922 x3=nextNode->pathlength()/2.;
00923
00924 if(x3==-1.) continue;
00925
00926
00927
00928
00929 if (x2> (x1+x3)) x2 = x2 - x1 - x3;
00930 else x2=0.;
00931
00932 cout
00933 <<"getTrackRadLength:"
00934 <<"\n\tIn Detector: "<<thisNode->getDetector()->getName()
00935 <<"\n\t\tMaterial: "<<thisNode->getDetector()->getMaterial()
00936 <<"\n\t\tLength: "<<x1
00937 <<"\t\tGap Length: "<<x2
00938 <<"\n\tNext Detector: "<<nextNode->getDetector()->getName()
00939 <<"\n\t\tMaterial: "<<nextNode->getDetector()->getMaterial()
00940 <<"\n\t\tLength: "<<x3
00941 << endl;
00942 {
00943 if (thisNode->getX0()>0) totalR += x1/thisNode->getX0();
00944 if (nextNode->getGasX0()>0) totalR += x2/nextNode->getGasX0();
00945 if (nextNode->getX0()>0) totalR += x3/nextNode->getX0();
00946 }
00947
00948 thisNode = nextNode;
00949 x1 = x3;
00950 }
00951 if (totalR>200.)
00952 cout <<"StiKalmanTrack::getTrackRadLength() -W- Total Rad Length Error: "<<totalR;
00953 return totalR;
00954 }
00955
00956 double StiKalmanTrack::getNearBeam(StThreeVectorD *pnt,StThreeVectorD *dir) const
00957 {
00958 StiKalmanTrackNode * inNode = lastNode;
00959 StThreeVectorD in(inNode->x_g(),inNode->y_g(),inNode->z_g());
00960
00961 StPhysicalHelixD hlx(fabs(inNode->getCurvature()),
00962 inNode->getDipAngle(),
00963 inNode->getPhase(),
00964 in,
00965 inNode->getHelicity());
00966 double per = hlx.period();
00967 double len = hlx.pathLength(0.,0.);
00968
00969
00970 if (fabs(len) > fabs(len+per)) len+=per;
00971 if (fabs(len) > fabs(len-per)) len-=per;
00972
00973 hlx.moveOrigin(len);
00974 if (pnt) (*pnt) = hlx.at(0);
00975
00976 if (dir) {
00977 double phase = hlx.phase();
00978 double dip = hlx.dipAngle();
00979 int h = hlx.h();
00980
00981 (*dir)[0]= -sin(phase)*cos(dip)*h;
00982 (*dir)[1]= cos(phase)*cos(dip)*h;
00983 (*dir)[2]= sin(dip)*h;}
00984
00985 return fabs(len);
00986 }
00987
00988
00999
01000 StiKalmanTrackNode * StiKalmanTrack::getInnOutMostNode(int inot,int qua) const
01001 {
01002 if (firstNode==0 || lastNode==0)
01003 {
01004
01005 throw runtime_error("StiKalmanTrack::getInnOutMostNode() -E- firstNode||lastNode==0");
01006 }
01007
01008 StiKalmanTrackNode *node;
01009 StiKTNIterator it =(inot) ? begin():rbegin();
01010 for (;(node=it());it++){
01011 if (!node->isValid()) continue;
01012 StiHit *hit = node->getHit();
01013 if ((qua&1) && !hit) continue;
01014 if ((qua&2) && hit && node->getChi2()>10000.) continue;
01015 return node;
01016 }
01017 cout << "StiKalmanTrack::getInnOutMostNode() -E- No requested nodes " << endl;
01018
01019 return 0;
01020 }
01021
01022 StiKalmanTrackNode * StiKalmanTrack::getOuterMostHitNode(int qua) const
01023 {
01024 return getInnOutMostNode(1,qua|1);
01025 }
01026
01027
01028
01039 StiKalmanTrackNode * StiKalmanTrack::getInnerMostHitNode(int qua) const
01040 {
01041 return getInnOutMostNode(0,qua|1);
01042 }
01043 #ifdef DO_TPCCATRACKER
01044 StiKalmanTrackNode * StiKalmanTrack::getInnerMostTPCHitNode(int qua) const
01045 {
01046 if (firstNode==0 || lastNode==0)
01047 {
01048
01049 throw runtime_error("StiKalmanTrack::getInnOutMostNode() -E- firstNode||lastNode==0");
01050 }
01051
01052 StiKalmanTrackNode *node = 0;
01053 StiKalmanTrackNode* leaf = getLastNode();
01054 StiKTNForwardIterator it(leaf);
01055 StiKTNForwardIterator end = it.end();
01056 for (;it!=end;++it)
01057 {
01058 StiKalmanTrackNode& node_t = *it;
01059 if (!node_t.isValid()) continue;
01060 if (node_t.getChi2()>10000.) continue;
01061 StiHit* hit = node_t.getHit();
01062 if (!hit) continue;
01063 if(hit->x()<58.f) continue;
01064 node = &node_t;
01065 return node;
01066 }
01067
01068 cout << "StiKalmanTrack::getInnOutMostNode() -E- No requested nodes " << endl;
01069
01070 return 0;
01071 }
01072 #endif
01073
01074 int StiKalmanTrack::getNNodes(int qua) const
01075 {
01076 StiKalmanTrackNode *node;
01077 StiKTNIterator it = begin();
01078 int nn=0;
01079 for (;(node=it());it++){
01080 if (!node->isValid()) continue;
01081 StiHit *hit = node->getHit();
01082 if ((qua&1) && !hit) continue;
01083 if ((qua&2) && hit && node->getChi2()>10000.) continue;
01084 nn++;
01085 }
01086 return nn;
01087 }
01088
01089
01091 vector<StiKalmanTrackNode*> StiKalmanTrack::getNodes(int detectorId) const
01092 {
01093 StiKTNIterator it;
01094 vector<StiKalmanTrackNode*> nodeVec;
01095 for (it=begin();it!=end();++it) {
01096 StiKalmanTrackNode* node = &(*it);
01097 const StiHit* hit = node->getHit();
01098 if(!hit) continue;
01099 const StiDetector *det = hit->detector();
01100 if (!det) continue;
01101 if (node->getDedx()<=0.) continue;
01102 if (detectorId!=det->getGroupId()) continue;
01103 nodeVec.push_back(node);
01104 }
01105 return nodeVec;
01106 }
01107
01108
01110 vector<const StMeasuredPoint*> StiKalmanTrack::stHits() const
01111 {
01112 StiKTNIterator it;
01113 vector<const StMeasuredPoint*> hits;
01114 for (it=begin();it!=end();++it) {
01115 const StiKalmanTrackNode* node = &(*it);
01116 if (!node->isValid()) continue;
01117 if (node->getChi2()>10000.) continue;
01118 const StiHit* hit = node->getHit();
01119 if (!hit) continue;
01120 if (!hit->detector()) continue;
01121 const StMeasuredPoint *stHit = hit->stHit();
01122 if (!stHit) continue;
01123 hits.push_back(stHit);
01124 }
01125 return hits;
01126 }
01127
01128
01129
01130
01138
01139 void StiKalmanTrack::reserveHits()
01140 {
01141 StiKTNForwardIterator it(lastNode);
01142 for_each( it, it.end(), SetHitUsed() );
01143 }
01144
01167
01168 StiTrackNode *StiKalmanTrack::extendToVertex(StiHit* vertex)
01169 {
01170 static int nCall=0; nCall++;
01171 double chi2;
01172 int dcaHit = vertex->isDca();
01173 StiKalmanTrackNode * sNode=0;
01174 StiKalmanTrackNode * tNode=0;
01175 bool trackExtended = false;
01176
01177 StiKalmanTrackNode * innerMostHitNode = getInnerMostHitNode();
01178 if (!innerMostHitNode) return 0;
01179
01180 StiHit localVertex = *vertex;
01181 sNode = getInnerMostNode();
01182 if (sNode->isDca()) {
01183 removeLastNode();
01184 sNode = getInnerMostNode();
01185 }
01186
01187
01188 localVertex.rotate(sNode->getAlpha());
01189 tNode = trackNodeFactory->getInstance();
01190 StiHit *myHit;
01191
01192
01193
01194
01195 if (tNode->propagate(sNode, &localVertex,kOutsideIn))
01196 {
01197
01198 double dy=0,dz=0,d=0;
01199 if (dcaHit) {
01200 tNode->setChi2(0);
01201 tNode->setHit(0);
01202 tNode->setDetector(0);
01203 return tNode;
01204 } else {
01205 tNode->setChi2(3e33);
01206 chi2 = tNode->evaluateChi2(&localVertex);
01207 dy=tNode->getY()- localVertex.y();
01208 dz=tNode->getZ()- localVertex.z();
01209 d = ::sqrt(dy*dy+dz*dz);
01210 _vChi2= chi2; _dca = d;
01211 }
01212
01213 #ifdef Sti_DEBUG
01214 int npoints[2] = {0,0};
01215 vector<StMeasuredPoint*> hitVec = stHits();
01216 for (vector<StMeasuredPoint*>::iterator point = hitVec.begin(); point!=hitVec.end();++point) {
01217 StHit * hit = dynamic_cast<StHit *>(*point);
01218 if (hit) {
01219 StDetectorId detId = hit->detector();
01220 if (detId == kTpcId) ++npoints[0];
01221 if (detId == kSvtId) ++npoints[1];
01222 }
01223 }
01224 cout << "StiKalmanTrack::extendToVertex: localVertex: " << localVertex << endl;
01225 cout << "StiKalmanTrack::extendToVertex: chi2 @ vtx: " << chi2
01226 << " dx:"<< dx
01227 << " dy:"<< dy
01228 << " dz:"<< dz
01229 << " d: "<< d
01230 << " dca: " << _dca << " npoints tpc/svt: " << npoints[0] << "/" << npoints[1] << endl;
01231 cout << "StiKalmanTrack::extendToVertex: TrackBefore:" << *this << endl;
01232 #endif
01233
01234
01235
01236 if (chi2<StiKalmanTrackFinderParameters::instance()->maxChi2Vertex())
01237 {
01238 myHit = StiToolkit::instance()->getHitFactory()->getInstance();
01239 *myHit = localVertex;
01240 tNode->setHit(myHit);
01241 tNode->setChi2(chi2);
01242 tNode->setDetector(0);
01243 trackExtended = (tNode->updateNode()==0);
01244
01245 #ifdef Sti_DEBUG
01246 cout << "StiKalmanTrack::extendToVertex: TrackAfter:" << *this << endl;
01247 #endif
01248 if (debug()) cout << "extendToVertex:: " << StiKalmanTrackNode::Comment() << endl;
01249
01250 if (trackExtended) return tNode;
01251 trackNodeFactory->free(tNode);
01252 }
01253 else if (d < 4) {
01254 LOG_DEBUG <<
01255 Form("Primary(%d) not accepted BUT d = %g chi2 = %g",nCall,d,chi2)
01256 << endm;
01257 }
01258 }
01259
01260
01261 return 0;
01262 }
01263
01264
01265 bool StiKalmanTrack::find(int direction)
01266 {
01267 static int nCall=0; nCall++;
01268
01269 StiDebug::Break(nCall);
01270 bool trackExtended=false;
01271 bool trackExtendedOut=false;
01272 int status = 0;
01273 setFlag(0);
01274
01275
01276 try
01277 {
01278 if (getNNodes(3)<4) return false;
01279 status = fit(kOutsideIn);
01280 if (getNNodes(3)<4) return false;
01281 if (debug()) cout << "StiKalmanTrack::find seed " << *((StiTrack *) this);
01282 double radSvt= 50.;
01283 if (trackFinder->find(this,kOutsideIn,radSvt)) {
01284 status = approx(1);
01285 status = refit() ; if(status) return false;
01286 trackExtended = getNNodes(3)>5;
01287 }
01288
01289 if (trackFinder->find(this,kOutsideIn,0.)) {
01290
01291 status = refit() ; if(status) return false;
01292 trackExtended = trackExtended || getNNodes(3)>5;
01293 }
01294 }
01295 catch (runtime_error & error)
01296 {
01297 cout << "SKT:find(int dir) -W- ERROR:" << error.what()<<endl;
01298 }
01299
01300 const StiKalmanTrackNode * outerMostNode = getOuterMostNode(2);
01301 if (!outerMostNode)
01302 {
01303 setFlag(-1);
01304 return false;
01305 }
01306 if (outerMostNode->getX()<185. )
01307 {
01308 try
01309 {
01310 if (debug()) cout << "StiKalmanTrack::find swap " << *((StiTrack *) this);
01311 trackExtendedOut= trackFinder->find(this,kInsideOut);
01312 if (trackExtendedOut) {
01313
01314 status = refit() ; if(status) return false;
01315 trackExtendedOut = getNNodes(3)>5;
01316 }
01317 if (debug()) cout << "StiKalmanTrack::find find(this,kInsideOut)" << *((StiTrack *) this);
01318 }
01319 catch (...)
01320 {
01321 cout << "StiKalmanTrack::find(int direction) -W- Exception while in insideOut find"<<endl;
01322 }
01323
01324 }
01325 setFlag(1);
01326
01327 return trackExtended||trackExtendedOut;
01328 }
01331
01332 vector<StiHit*> StiKalmanTrack::getHits()
01333 {
01334 vector<StiHit*> hits;
01335 StiKalmanTrackNode* leaf = getLastNode();
01336 StiKTNForwardIterator it(leaf);
01337 StiKTNForwardIterator end = it.end();
01338 for (;it!=end;++it)
01339 {
01340 const StiKalmanTrackNode& node = *it;
01341 if (!node.isValid()) continue;
01342 if (node.getChi2()>10000.) continue;
01343 StiHit* hit = node.getHit();
01344 if (!hit) continue;
01345 hits.push_back(hit);
01346 }
01347 return hits;
01348 }
01349
01350
01352 double StiKalmanTrack::getDca(const StiHit * vertex) const
01353 {
01354 StiKalmanTrackNode* node;
01355
01356 node = getInnerMostHitNode(2);
01357 StThreeVectorD originD(node->x_g(),node->y_g(),node->z_g());
01358 StThreeVectorD vxDD(vertex->x_g(), vertex->y_g(),vertex->z_g());
01359 StPhysicalHelixD physicalHelix(0.,0.,0.,originD,-1);
01360 physicalHelix.setParameters(fabs(node->getCurvature()),
01361 node->getDipAngle(),
01362 node->getPhase(),
01363 originD,
01364 node->getHelicity());
01365 double dca = physicalHelix.distance(vxDD);
01366 return dca;
01367 }
01368
01369 ostream& operator<<(ostream& os, const StiKalmanTrack& track)
01370 {
01371 try
01372 {
01373 os << *((StiTrack *) &track);
01374 os <<"List of nodes" << endl;
01375 StiKTNIterator tNode = track.begin();
01376 StiKTNIterator eNode = track.end();
01377
01378
01379 while (tNode != eNode) {
01380 StiKalmanTrackNode *thisNode = &(*tNode);
01381 if (thisNode) os << *thisNode;
01382 tNode++;
01383 }
01384 }
01385 catch (runtime_error & rte)
01386 {
01387 os << " Run-time Error while accessing track parameters: " << rte.what() << endl;
01388 }
01389 catch (logic_error & le)
01390 {
01391 os << " Logic Error while accessing track parameters: " << le.what() << endl;
01392 }
01393 return os;
01394 }
01395
01396
01401 StiKalmanTrackNode * StiKalmanTrack::extrapolateToBeam()
01402 {
01403 StiKalmanTrackNode * innerMostNode = getInnerMostNode();
01404
01405 if (!innerMostNode) return 0;
01406 StiKalmanTrackNode * n = trackNodeFactory->getInstance();
01407 if (n->propagateToBeam(innerMostNode,kOutsideIn)) return n;
01408 trackNodeFactory->free(n);
01409 return 0;
01410 }
01411
01412
01413 StiKalmanTrackNode * StiKalmanTrack::extrapolateToRadius(double radius)
01414 {
01415 StiKalmanTrackNode * outerMostNode = getOuterMostNode();
01416
01417 if (!outerMostNode) return 0;
01418 StiKalmanTrackNode *n = trackNodeFactory->getInstance();
01419 if (n->propagateToRadius(outerMostNode,radius,kOutsideIn)) return n;
01420 trackNodeFactory->free(n);
01421 return 0;
01422 }
01423
01424
01425 void StiKalmanTrack::add(StiTrackNode * node,int direction)
01426 {
01427
01428 StiKalmanTrackNode *Node = (StiKalmanTrackNode*)node;
01429 if (lastNode==0) {
01430 lastNode = firstNode = Node; return;
01431 }
01432 if (direction==0) {
01433 lastNode->add(Node,direction);
01434 lastNode = Node;
01435 } else {
01436 firstNode->add(Node,direction);
01437 firstNode = Node;
01438 }
01439 }
01440
01441 void StiKalmanTrack::setFirstLastNode(StiKalmanTrackNode * node)
01442 {
01443 firstNode = (StiKalmanTrackNode*)node->getFirstNode();
01444 lastNode = (StiKalmanTrackNode*)node->getLastNode ();
01445 }
01446
01447 void StiKalmanTrack::removeLastNode()
01448 {
01449 StiKalmanTrackNode *node = lastNode;
01450 lastNode = (StiKalmanTrackNode*)node->disconnect();
01451 BFactory::Free(node);
01452 }
01453
01454 int StiKalmanTrack::refit()
01455 {
01456 int errType = kNoErrors;
01457
01458 static int nCall=0; nCall++;
01459 StiDebug::Break(nCall);
01460 enum {kMaxIter=30,kPctLoss=10,kHitLoss=3};
01461 static double defConfidence = StiDebug::dFlag("StiConfidence",0.01);
01462 int nNBeg = getNNodes(3), nNEnd = nNBeg;
01463 if (nNBeg<=3) return 1;
01464 if (!mgMaxRefiter) return 0;
01465 StiKalmanTrackNode *inn= getInnerMostNode(3);
01466 int fail=0,status=0;
01467
01468 StiNodePars pPrev;
01469 StiNodeErrs ePrev;
01470 int iter=0,igor=0;
01471 double qA;
01472 double errConfidence = defConfidence;
01473 for (int ITER=0;ITER<mgMaxRefiter;ITER++) {
01474 for (iter=0;iter<kMaxIter;iter++) {
01475 fail = 0;
01476 errType = kNoErrors;
01477 sTNH.set(StiKalmanTrackFitterParameters::instance()->getMaxChi2()*10,StiKalmanTrackFitterParameters::instance()->getMaxChi2Vtx()*100,errConfidence,iter);
01478 pPrev = inn->fitPars();
01479 ePrev = inn->fitErrs();
01480
01481 status = refitL();
01482 if (status) {fail= 1; errType = kRefitFail; break;}
01483 nNEnd = sTNH.getUsed();
01484 if ((nNEnd <=3)) {fail= 2; errType = kNotEnoughUsed; break;}
01485 if (!inn->isValid() || inn->getChi2()>1000) {
01486 inn = getInnerMostNode(3); fail=-1; errType = kInNodeNotValid; continue;}
01487 qA = diff(pPrev,ePrev,inn->fitPars(),inn->fitErrs(),igor);
01488 static int oldRefit = StiDebug::iFlag("StiOldRefit");
01489 if (oldRefit) {
01490 if (qA>0.5) {fail=-2; errType = kBadQA; continue;}
01491 } else {
01492 if (qA <1 && errConfidence>0.1) errConfidence = 0.1;
01493 if (qA>0.01) {fail=-2; errType = kBadQA; continue;}
01494 if (sTNH.isCutStep()) {fail=-2; errType = kBadQA; continue;}
01495 }
01496 double info[2][8];
01497 sTNH.mCurvQa.getInfo(info[0]);
01498 sTNH.mTanlQa.getInfo(info[1]);
01499 break;
01500 }
01501 if (fail>0) break;
01502
01503 StiKalmanTrackNode *worstNode= sTNH.getWorst();
01504 if (worstNode && worstNode->getChi2()>StiKalmanTrackFitterParameters::instance()->getMaxChi2())
01505 {
01506 worstNode->setHit(0); worstNode->setChi2(3e33); continue;}
01507 if (rejectByHitSet()) { releaseHits() ; continue;}
01508
01509 if (!fail) break;
01510
01511 StiKalmanTrackNode *flipFlopNode= sTNH.getFlipFlop();
01512 if (flipFlopNode && flipFlopNode->getFlipFlop()>kMaxIter/3)
01513 {
01514 flipFlopNode->setHit(0); flipFlopNode->setChi2(3e33); continue;}
01515 break;
01516
01517
01518
01519 }
01520 StiKalmanTrackNode *vertexNode= sTNH.getVertexNode();
01521
01522
01523 while (!fail && vertexNode) {
01524 fail = 13;
01525 errType = kVertexNodeInvalid;
01526 if (!vertexNode->isValid()) break;
01527 fail = 99;
01528 errType = kNodeNotValid;
01529 if ( vertexNode->getChi2()>StiKalmanTrackFitterParameters::instance()->getMaxChi2Vtx()) break;
01530 fail = 98;
01531 errType = kTooManyDroppedNodes;
01532 if (nNBeg*kPctLoss/100 < nNBeg-nNEnd
01533 && nNEnd+kHitLoss < nNBeg) break;
01534 fail = 0;
01535 errType = kNoErrors;
01536 break;
01537 }
01538 if (!fail) {
01539 StiKalmanTrackNode *node;
01540 StiKTNIterator it = begin();
01541 for (;(node=it());it++){
01542 if (node == vertexNode) continue;
01543 StiHit *hit = node->getHit();
01544 if(!hit) continue;
01545 hit->setTimesUsed(0);
01546 node->setHit(0);
01547 if (!node->isValid()) continue;
01548 if (node->getChi2()>10000.) continue;
01549 assert(node->getChi2()<=StiKalmanTrackFitterParameters::instance()->getMaxChi2());
01550 hit->setTimesUsed(1);
01551 node->setHit(hit);
01552 }
01553 }
01554 static int VPDEBUG=0;
01555 if (VPDEBUG) {
01556 if (fail>0) {
01557 LOG_DEBUG <<
01558 Form("StiKalmanTrack::refit(%d)=%d ***FAILED*** %d>%d",fail,nCall,nNBeg,nNEnd)
01559 << endm;
01560 } else {
01561 LOG_DEBUG <<
01562 Form("StiKalmanTrack::refit(%d)=%d ***EndIters*** %d>%d",fail,nCall,nNBeg,nNEnd)
01563 << endm;
01564 }
01565 }
01566
01567 if (fail) setFlag(-1);
01568 #ifdef DO_TPCCATRACKER
01569 return errType;
01570 #else
01571 return fail;
01572 #endif
01573 }
01574
01575 int StiKalmanTrack::refitL()
01576 {
01577 static int nCall=0;nCall++;
01578 StiDebug::Break(nCall);
01579
01580 StiKTNIterator source;
01581 StiKalmanTrackNode *pNode = 0,*targetNode;
01582 int iNode=0, status = 0,isStarted=0,restIsWrong=0;
01583 for (source=rbegin();source!=rend();source++) {
01584 iNode++;
01585 targetNode = &(*source);
01586 if (restIsWrong) { targetNode->setInvalid(); continue;}
01587
01588 if (!isStarted) {
01589 if (!targetNode->getHit()) targetNode->setInvalid();
01590 if ( targetNode->getChi2()>1000) targetNode->setInvalid();
01591 if (!targetNode->isValid()) continue;
01592 }
01593 isStarted++;
01594 sTNH.set(pNode,targetNode);
01595 status = sTNH.makeFit(0);
01596 if (status) {restIsWrong = 2005; targetNode->setInvalid();}
01597 if (!targetNode->isValid()) continue;
01598 pNode = targetNode;
01599 }
01600
01601 pNode = 0; iNode=0;isStarted=0;restIsWrong=0;
01602 for (source=begin();source!=end();source++) {
01603 iNode++;
01604 targetNode = &(*source);
01605 if (restIsWrong) { targetNode->setInvalid(); continue;}
01606 if (!isStarted) {
01607 if (!targetNode->getHit()) targetNode->setInvalid();
01608 if ( targetNode->getChi2()>1000) targetNode->setInvalid();
01609 if (!targetNode->isValid()) continue;
01610 }
01611 isStarted++;
01612 sTNH.set(pNode,targetNode);
01613 status = sTNH.makeFit(1);
01614 if (status) {restIsWrong = 2005; targetNode->setInvalid();}
01615 if (!targetNode->isValid()) continue;
01616 pNode = targetNode;
01617 }
01618 return 0;
01619 }
01620
01621 void StiKalmanTrack::reduce()
01622 {
01623 StiKTNIterator source;
01624 for (source=begin();source!=end();source++) {(*source).reduce();}
01625 }
01626
01627 void StiKalmanTrack::unset()
01628 {
01629 if (!lastNode) return;
01630 StiKTNIterator source;
01631 for (source=begin();source!=end();source++) {BFactory::Free(&(*source));}
01632 lastNode=0; firstNode=0;
01633 }
01634
01635 void StiKalmanTrack::print(const char *opt) const
01636 {
01637 printf("Track %p",(void*)this);
01638
01639 StiKTNIterator it;
01640 int n=0;
01641 for (it=begin();it!=end();++it) {
01642 StiKalmanTrackNode *node = &(*it);
01643 StiHit *hit = node->getHit();
01644 if (!hit && strchr(opt,'h')) continue;
01645 if (!hit && strchr(opt,'H')) continue;
01646 n++;
01647 printf("%3d - ",n);
01648 node->print(opt);
01649 }
01650 }
01651
01652 #ifdef APPROX_DEBUG
01653 #include "TCanvas.h"
01654 #include "TH1F.h"
01655 #include "TProfile.h"
01656 #endif // APPROX_DEBUG
01657
01658 int StiKalmanTrack::approx(int mode)
01659 {
01660 static int nCall=0; nCall++;
01661 StiDebug::Break(nCall);
01662
01663 #ifdef APPROX_DEBUG
01664 static TCanvas *myCanvas=0;
01665 static TH1 *H[4];
01666
01667 if(!myCanvas) {
01668 myCanvas=new TCanvas("Approx","",600,800);
01669 H[0] = new TH1F("Approx0l","Approx0l", 100,0,5);
01670 H[1] = new TH1F("Approx1l","Approx1l", 100,0,5);
01671 H[2] = new TProfile("Approx0 ","Approx0", 30,0,30);
01672 H[3] = new TProfile("Approx1 ","Approx1 ", 30,0,30);
01673 myCanvas->Divide(1,4);
01674 for (int i=0;i<4;i++) {myCanvas->cd(i+1); H[i]->Draw();}
01675 }
01676 #endif // APPROX_DEBUG
01677
01678 const double BAD_XI2[2]={99,22},XI2_FACT=9;
01679 int nNode,nNodeIn,iNode=0;
01680 double Xi2=0;
01681 StiHitErrs hr;
01682
01683
01684 StiKTNIterator source;
01685 StiKalmanTrackNode *targetNode;
01686 nNode=0;
01687 double hz=0;
01688 THelixFitter circ;
01689 THelixTrack cirl;
01690 for (source=rbegin();(targetNode=source());++source) {
01691 iNode++;
01692 if (!hz) hz = targetNode->getHz();
01693 if (!targetNode->isValid()) continue;
01694 const StiHit * hit = targetNode->getHit();
01695 if (!hit) continue;
01696 if (targetNode->getChi2()>1000) continue;
01697 circ.Add(hit->x_g(),hit->y_g(),hit->z_g());
01698 hr = targetNode->getGlobalHitErrs(hit);
01699 circ.AddErr(hr.A,hr.hZZ);
01700 nNode++;
01701 }
01702 if (!nNode) return 1;
01703 nNodeIn = nNode;
01704 int nPnts = nNode;
01705 if (nPnts==2) {
01706 nPnts=3;
01707 circ.Add(0.,0.,0.);
01708 double vErr[3]={1.,0.,1};
01709 circ.AddErr(vErr,100*100.);
01710 }
01711
01712
01713 Xi2 =circ.Fit();
01714 if (mode==1 && Xi2>BAD_XI2[1]) return 2;
01715 #ifdef APPROX_DEBUG
01716 H[mode+0]->Fill(log(Xi2)/log(10.));
01717 H[mode+2]->Fill(nNode,Xi2);
01718 #endif // APPROX_DEBUG
01719
01720 circ.MakeErrs();
01721
01722 double s=0,xyz[3];
01723 double curv = circ.GetRho();
01724 iNode = 0;
01725 for (source=rbegin();(targetNode=source());++source) {
01726 iNode++;
01727 if (!targetNode->isValid()) continue;
01728 const StiHit *hit = targetNode->getHit();
01729 if (hit) {
01730 xyz[0] = hit->x_g();
01731 xyz[1] = hit->y_g();
01732 xyz[2] = hit->z_g();
01733 } else {
01734 xyz[0] = targetNode->x_g();
01735 xyz[1] = targetNode->y_g();
01736 xyz[2] = targetNode->z_g();
01737 }
01738 double ds = circ.Path(xyz[0],xyz[1]);
01739 circ.Move(ds);
01740 s+=ds;
01741 cirl = circ;
01742 double alfa = targetNode->getAlpha();
01743 cirl.Rot(-alfa);
01744 StiNodePars P = targetNode->fitPars();
01745 P.x() = cirl.Pos()[0];
01746 P.y() = cirl.Pos()[1];
01747 P.z() = cirl.Pos()[2];
01748 P.eta() = atan2(cirl.Dir()[1],cirl.Dir()[0]);
01749 P.curv() = curv;
01750 double hh = P.hz();
01751 hh = (fabs(hh)<1e-10)? 0:1./hh;
01752 P.ptin() = (hh)? curv*hh:1e-3;
01753 P.tanl() = cirl.GetSin()/cirl.GetCos();
01754 P._cosCA = cirl.Dir()[0]/cirl.GetCos();
01755 P._sinCA = cirl.Dir()[1]/cirl.GetCos();
01756 targetNode->fitPars() = P;
01757 int ians = targetNode->nudge();
01758 if(ians) {nNode--; targetNode->setInvalid();continue;}
01759 P = targetNode->fitPars();
01760 StiNodeErrs &E = targetNode->fitErrs();
01761 cirl.StiEmx(E.A);
01762 TCL::vscale(&(E._cPX),hh,&(E._cPX),5);
01763 E._cPP*=hh; E._cTP*=hh;
01764 if ((mode&1)==0 && Xi2>XI2_FACT) E*=Xi2/XI2_FACT;
01765 E.check("In aprox");
01766 }
01767
01768 if (Xi2>BAD_XI2[mode])return 2;
01769 if (nNode==nNodeIn) return 0;
01770 if (nNode<2) return 3;
01771 return 0;
01772 }
01773
01774 double diff(const StiNodePars &p1,const StiNodeErrs &e1
01775 ,const StiNodePars &p2,const StiNodeErrs &e2,int &igor)
01776 {
01777 double est=0;
01778 for (int i=0;i<kNPars;i++) {
01779 double err = 0.5*(e1(i,i)+e2(i,i));
01780 if (err<1e-10) continue;
01781 double dif = pow(p1[i]-p2[i],2)/err;
01782 if (est<dif) {est = dif; igor = i;}
01783 }
01784 return est;
01785 }
01786
01787 StiKalmanTrack &StiKalmanTrack::operator=(const StiKalmanTrack &tk)
01788 {
01789 StiTrack::operator=(tk);
01790 firstNode=0;
01791 lastNode=0;
01792
01793 mSeedHitCount=tk.mSeedHitCount;
01794 mFlag =tk.mFlag;
01795 m =tk.m;
01796 _dca =tk._dca;
01797 _vChi2 =tk._vChi2;
01798 mVertex =tk.mVertex;
01799 StiKTNIterator it;
01800 for (it=tk.begin();it!=tk.end();it++){
01801 const StiKalmanTrackNode *node = &(*it);
01802 if (!node->isValid()) continue;
01803 StiKalmanTrackNode *myNode=trackNodeFactory->getInstance();
01804 *myNode=*node;
01805 add(myNode,kOutsideIn);
01806 }
01807 return *this;
01808 }
01809
01810
01811 StThreeVector<double> StiKalmanTrack::getPoint(int firstLast) const
01812 {
01813 const StiKalmanTrackNode* node = getInnOutMostNode(firstLast,3);
01814 return StThreeVector<double>(node->x_g(),node->y_g(),node->z_g());
01815 }
01816
01817
01818 void StiKalmanTrack::setMaxRefiter(int maxRefiter)
01819 {
01820 mgMaxRefiter = maxRefiter;
01821 }
01822
01823 int StiKalmanTrack::rejectByHitSet() const
01824 {
01825 StiKalmanTrackNode *node;
01826 int sum=0;
01827 for (StiKTNIterator it = rbegin();(node=it());it++){
01828 if (node->x()>50) break;
01829 if (!node->isValid()) continue;
01830 StiHit *hit = node->getHit();
01831 if (!hit) continue;
01832 if (!hit->detector()) continue;
01833 if (node->getChi2()>1000) continue;
01834 sum+= StiKalmanTrackFinderParameters::instance()->hitWeight(int(hit->x()));
01835 }
01836 if (!sum) return 0;
01837 return sum < StiKalmanTrackFinderParameters::instance()->sumWeight();
01838 }
01839
01840 int StiKalmanTrack::releaseHits(double rMin,double rMax)
01841 {
01842 StiKalmanTrackNode *node;
01843 int sum=0;
01844 for (StiKTNIterator it = rbegin();(node=it());it++){
01845 StiHit *hit = node->getHit();
01846 if (!hit) continue;
01847 if (!hit->detector()) continue;
01848 if (hit->x()<rMin) continue;
01849 if (hit->x()>rMax) break;
01850 sum++;
01851 node->setHit(0);
01852 hit->setTimesUsed(0);
01853 }
01854 return sum;
01855 }
01856
