#include <MdcTrackListCsmc.h>
Inheritance diagram for MdcTrackListCsmc:
Public Member Functions | |
int | arbitrateHits () |
int | arbitrateHits () |
int | createFromSegs (MdcSegList *segs, const MdcHitMap *, const MdcDetector *, TrkContext &, double trackT0) |
int | createFromSegs (MdcSegList *segs, const MdcHitMap *, const MdcDetector *, TrkContext &, double trackT0) |
void | dropMultiHotInLayer (const MdcTrack *tk) |
void | dropMultiHotInLayer (const MdcTrack *tk) |
int | finish3d (TrkRecoTrk &trk) |
int | finish3d (TrkRecoTrk &trk) |
MdcTrackListCsmc (const MdcTrackParams &tkPar) | |
MdcTrackListCsmc (const MdcTrackParams &tkPar) | |
void | newParams (const MdcTrackParams &tkPar) |
void | newParams (const MdcTrackParams &tkPar) |
int | nTrack () const |
int | nTrack () const |
void | plot () const |
void | plot () const |
void | remove (MdcTrack *atrack) |
void | remove (MdcTrack *atrack) |
void | setD0Cut (double d0Cut) |
void | setD0Cut (double d0Cut) |
void | setPlot (int plotFlag) |
void | setPlot (int plotFlag) |
void | setPtCut (double ptCut) |
void | setPtCut (double ptCut) |
void | setZ0Cut (double z0Cut) |
void | setZ0Cut (double z0Cut) |
void | store (RecMdcTrackCol *, RecMdcHitCol *) |
void | store (RecMdcTrackCol *, RecMdcHitCol *) |
~MdcTrackListCsmc () | |
~MdcTrackListCsmc () | |
Static Public Attributes | |
double | m_d0Cut = -999. |
double | m_ptCut = -999. |
double | m_z0Cut = -999. |
Protected Attributes | |
MdcTrackParams | tkParam |
|
00034 : 00035 MdcTrackListBase(tkPar) { 00036 // ************************************************************************* 00037 return; 00038 }
|
|
00041 {}
|
|
|
|
|
|
|
|
00104 { 00105 //************************************************************************* 00106 // Look at all hits used in two or more tracks. Assign hits to the track 00107 // that gives the lower residual. If, however, many hits are shared by 00108 // a pair of tracks, assign them all to one or the other. 00109 // Refit any tracks that have had hits dropped. 00110 // The implementation is very clumsy, since the arrays were originally 00111 // indexed by id # => there is an unneeded layer of indexing. 00112 00113 // return # of tracks deleted 00114 00115 if (8 == tkParam.lPrint){ 00116 std::cout << "=======Print before arbitrateHits=======" << std::endl; 00117 } 00118 00119 int nDeleted = 0; 00120 std::vector<MdcTrack*> trksToKill; 00121 trksToKill.reserve(4); 00122 00123 MdcMap<long,long> idMap; 00124 00125 //usedInTrackNum records how many shared hits track has with each other track 00126 int* usedInTrackNum = new int [nTrack()]; 00127 // to navigate from track id # to track pointer: 00128 MdcTrack** trkXRef = new MdcTrack* [nTrack()]; 00129 //refitTrack flags track id #s of tracks to be refit 00130 int *refitTrack = new int [nTrack()]; 00131 for (int i = 0; i < nTrack(); i++) { 00132 refitTrack[i] = 0; 00133 } 00134 00135 // Fill xref table 00136 int itrack; 00137 for (itrack = 0; itrack < nTrack(); itrack++) { 00138 MdcTrack *atrack = (*this)[itrack]; 00139 if (atrack == 0) continue; // I don't think it can be, but . . . 00140 idMap.put(atrack->track().id(), itrack); 00141 trkXRef[itrack] = atrack; 00142 } 00143 // Loop through the tracks 00144 for (itrack = 0; itrack < nTrack(); itrack++) { 00145 00146 if (8 == tkParam.lPrint) std::cout<<"arbitrate track No."<<itrack<< std::endl; 00147 MdcTrack *atrack = (*this)[itrack]; 00148 if (atrack == 0) continue; 00149 TrkRecoTrk& aRecoTrk = atrack->track(); 00150 int lRefit = 0; 00151 int trackOld = -1; 00152 const TrkFit* tkFit = aRecoTrk.fitResult(); 00153 assert (tkFit != 0); 00154 TrkHitList* hitList = aRecoTrk.hits(); 00155 assert (hitList != 0); 00156 restart: 00157 for (int ii = 0; ii < nTrack(); ii++) usedInTrackNum[ii] = 0; 00158 00159 // Loop through hits on track, counting # used in other tracks 00160 int nPrev = 0; 00161 int nHitDeleted = 0; 00162 int nHitInLayer[43];//yzhang 2010-09-20 for bad tracking testing 00163 int nDeleteInLayer[43];//yzhang 2010-09-20 00164 for(int i=0;i<43;i++){ 00165 nHitInLayer[i]=0; 00166 nDeleteInLayer[i]=0; 00167 } 00168 if(8 == tkParam.lPrint) std::cout<< "--arbitrate--"<<std::endl; 00169 for (TrkHitList::hot_iterator ihit(hitList->begin()); ihit != hitList->end(); ++ihit){ 00170 int nUsed = ihit->hit()->nUsedHits(); 00171 if (8 == tkParam.lPrint){ 00172 std::cout<<"nUsed="<<nUsed<<":"; 00173 ihit->hit()->printAll(std::cout); 00174 } 00175 if (8 == tkParam.lPrint) { 00176 double deltaChi = -999; 00177 ihit->getFitStuff(deltaChi); 00178 std::cout<< "deltaChi="<<deltaChi<<std::endl; 00179 } 00180 int layer = ihit->layerNumber(); 00181 nHitInLayer[layer]++; 00182 00183 if (!ihit->isActive()) { 00184 //----------------------------------- 00185 //yzhang delete not ACT hit 2010-05-14 00186 //----------------------------------- 00187 if(tkParam.lRemoveInActive ) {//2010-05-16 00188 if (8 == tkParam.lPrint) { 00189 std::cout<< "=remove above inactive "<<std::endl; 00190 } 00191 TrkFundHit* hit = const_cast<TrkFundHit*> (ihit->hit()); 00192 if(hitList->removeHit(hit)){ hit->setUnusedHit(&(*ihit)); } 00193 if(ihit == hitList->end()) break; 00194 --ihit;//be careful of the iterator, yzhang 00195 nDeleteInLayer[layer]++; 00196 } 00197 continue; // active hits only yzhang 2009-11-03 delete 00198 } 00199 if (nUsed > 1) { 00200 bool wasUsed = false; 00201 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator> q = 00202 ihit->hit()->getUsedHits(); 00203 for (TrkFundHit::hot_iterator i = q.first; i != q.second; ++i) { 00204 if ( !i->isActive() ) continue; //yzhang 2009-11-03 delete 00205 TrkRecoTrk * recoTrk=i->parentTrack(); 00206 int id = recoTrk->id(); 00207 if (id == aRecoTrk.id()) continue; //skip same track 00208 long index = 0; 00209 idMap.get(id, index); 00210 assert(index >= 0); 00211 usedInTrackNum[index]++; 00212 if (8 == tkParam.lPrint){ 00213 std::cout<<" track "<<itrack<<"&" <<index 00214 << " shared hits "<<usedInTrackNum[index]<<":"; 00215 ihit->printAll(std::cout); 00216 } 00217 wasUsed = true; 00218 } 00219 if (wasUsed) nPrev++; 00220 }// end nUsed > 1 00221 } // end loop over hits 00222 00223 for (int i=0;i<43;i++){ 00224 if (8 == tkParam.lPrint) { 00225 std::cout<<i<<" nHitInLayer "<<nHitInLayer[i] 00226 <<" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl; 00227 } 00228 if(nHitInLayer[i]>0 && (nHitInLayer[i]-nDeleteInLayer[i])==0) { 00229 nHitDeleted++; 00230 if (8 == tkParam.lPrint) { 00231 cout << "rec hits have been deleted in this layer"<<std::endl; 00232 } 00233 } 00234 } 00235 00236 if (8 == tkParam.lPrint) { cout << " nHitDeleted = "<<nHitDeleted<<std::endl; } 00237 //yzhang add nHitDeleted cut 2010-09-13 00238 // remove track if # not Active 00239 if (nHitDeleted >= tkParam.nHitDeleted) { 00240 if (8 == tkParam.lPrint) { 00241 cout << " nHitDeleted "<<nHitDeleted<<" > "<<tkParam.nHitDeleted 00242 <<" Killing tkNo " << itrack << endl; 00243 } 00244 nDeleted++; 00245 delete &(atrack->track()); // Delete the RecoTrk inside atrack 00246 atrack->setTrack(0); 00247 trksToKill.push_back(atrack); 00248 continue; 00249 } 00250 00251 //******* 00252 // How many hits are shared with a single track? 00253 int nMost = 0; 00254 int trackMost = 0; 00255 for (int ii = 0; ii < nTrack(); ii++) { 00256 if (8 == tkParam.lPrint){ 00257 std::cout<<"tk:"<<itrack<<"&"<<ii 00258 <<" shared "<<usedInTrackNum[ii]<<" hits "<< std::endl; 00259 } 00260 if (usedInTrackNum[ii] > nMost) { 00261 nMost = usedInTrackNum[ii]; 00262 trackMost = ii; //index of track w/ most hits in common w/ current trk 00263 } 00264 } 00265 00266 // A little precaution against infinite loops: 00267 if (trackMost == trackOld) { 00268 std::cout << "ErrMsg(error) MdcTrackListBase:" 00269 << "Something ghastly happened in MdcTrackListBase::arbitrateHits" 00270 << std::endl; 00271 return 0; 00272 } 00273 trackOld = trackMost; 00274 00275 00276 //****** 00277 // Decide whether to handle hits individually or in group 00278 double groupDiff = 0.0; // relative quality of grouped hits for the two 00279 // tracks; > 0. => current track worse 00280 int nFound = 0; // # of grouped hits located so far 00281 TrkHitOnTrk **theseHits = 0; // grouped hits as seen in current track 00282 TrkHitOnTrk **thoseHits = 0; // grouped hits as seen in the other track 00283 int lGroupHits = 0; 00284 00285 if (nMost >= tkParam.nOverlap) { 00286 if (8 == tkParam.lPrint){ 00287 std::cout<<"track "<<trackMost<<" shared "<<nMost<<" hits > Cut nOverlap " 00288 <<tkParam.nOverlap<<", group hits!"<<std::endl; 00289 } 00290 lGroupHits = 1; 00291 theseHits = new TrkHitOnTrk*[nMost]; 00292 thoseHits = new TrkHitOnTrk*[nMost]; 00293 } 00294 00295 //********* 00296 // Go back through hits on this track, looking up the overlap of each 00297 // if grouping hits, only deal with hits shared with trackMost on this pass 00298 // otherwise, deal with all shared hits as encountered 00299 if(8 == tkParam.lPrint) std::cout<<"Go back through hits, looking up overlap hits"<< std::endl; 00300 if (nMost > 0) { 00301 if (8 == tkParam.lPrint) std::cout<<" nHits= "<< hitList->nHit()<< std::endl; 00302 for (TrkHitList::hot_iterator ihit(hitList->begin()); ihit != hitList->end(); ++ihit) { 00303 int nUsed = ihit->hit()->nUsedHits(); 00304 00305 if (8 == tkParam.lPrint){ 00306 std::cout<< "--hit go back, nUsed="<<nUsed<<":"; 00307 ihit->hit()->printAll(std::cout); 00308 } 00309 00310 // only shared hits 00311 if (nUsed < 2) { continue; } 00312 00313 // active hits only 00314 if (!ihit->isActive()) { 00315 if (8 == tkParam.lPrint){ std::cout<<"act=0 continue"<<std::endl; } 00316 continue; 00317 } 00318 00319 //*** look at all overlaps for this hit 00320 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator> q = ihit->hit()->getUsedHits(); 00321 while (q.first!=q.second) { // nUsed > 0 00322 int dropThisHit = 0; 00323 TrkHitOnTrk *otherHot = const_cast<TrkHitOnTrk*>((--q.second).get()); 00324 TrkRecoTrk *otherTrack = otherHot->parentTrack(); 00325 00326 if (!otherHot->isActive()) continue; 00327 00328 // Again, skip "overlap" of track with itself 00329 if ( &aRecoTrk == otherTrack) continue; 00330 int otherId = otherTrack->id(); 00331 long otherIndex = -1; 00332 idMap.get(otherId, otherIndex); assert(otherIndex >= 0); 00333 00334 // if grouping hits, only look at hits shared with trackMost 00335 if (lGroupHits && otherIndex != trackMost) continue; 00336 00337 if (lGroupHits) { 00338 if (8 == tkParam.lPrint) { 00339 std::cout<<"group hits "<< std::endl; 00340 } 00341 // Calculate contribution of group to each chisq/dof 00342 // groupDiff += fabs(ihit->resid(0)) - 00343 // fabs(otherHot->resid(0)); 00344 // Hack to handle tracks with 5 active hits: 00345 int aDof = tkFit->nActive() - 5; 00346 assert (otherTrack->fitResult() != 0); 00347 int otherDof = otherTrack->fitResult()->nActive() - 5; 00348 if (aDof <= 0) {groupDiff = 999;} 00349 else if (otherDof <= 0) {groupDiff = -999;} 00350 else { 00351 groupDiff += ihit->resid(0) * ihit->resid(0) * ihit->weight() / 00352 aDof - 00353 otherHot->resid(0) * otherHot->resid(0) * otherHot->weight() / 00354 otherDof; 00355 } 00356 theseHits[nFound] = const_cast<TrkHitOnTrk*>(ihit.get()); 00357 thoseHits[nFound] = otherHot; 00358 nFound++; 00359 dropThisHit = 1; 00360 } else { // handle hits individually 00361 00362 if (8 == tkParam.lPrint) { 00363 std::cout<<"handle hits individually"<< std::endl; 00364 } 00365 nFound++; 00366 if (fabs(ihit->resid(0)) > fabs(otherHot->resid(0)) ) { 00367 // turn off (inactivate) hit on this track 00368 lRefit = 1; 00369 // ihit->hit()->setUnusedHit(ihit.get()); 00370 //Should I be setting inactive, or deleting the hit??????? 00371 const_cast<TrkHitOnTrk*>(ihit.get())->setActivity(0); 00372 dropThisHit = 1; 00373 if (8 == tkParam.lPrint) { 00374 std::cout<<"dorp hit "; 00375 const_cast<TrkHitOnTrk*>(ihit.get())->print(std::cout); 00376 } 00377 break; // found other hit, so quit loop 00378 } else { 00379 // inactivate hit on other track 00380 refitTrack[otherIndex] = 1; 00381 // otherHot->hit()->setUnusedHit(otherHot); 00382 otherHot->setActivity(0); 00383 if (8 == tkParam.lPrint) { 00384 std::cout<<"inactive hit on other track"; 00385 const_cast<TrkHitOnTrk*>(ihit.get())->print(std::cout); 00386 } 00387 break; // found other hit, so quit loop 00388 } 00389 } // end grouped/individual treatment 00390 00391 if (dropThisHit == 1) break; // don't look for other matches since 00392 // this hit is now turned off 00393 } // end loop over nUsed 00394 00395 // Quit if we've found all of the shared hits on this track 00396 if (lGroupHits && nFound == nMost || nFound == nPrev) { 00397 if (8 == tkParam.lPrint) { 00398 std::cout<<"we've found all of the shared hits on this track,Quit"<<std::endl; 00399 } 00400 break; 00401 } 00402 00403 } // end loop over hits 00404 00405 // Decide which track grouped hits belong with and inactivate accordingly 00406 if (lGroupHits) { 00407 if (8 == tkParam.lPrint) { 00408 cout << "nGroup: " << nMost << " groupDiff: " << groupDiff << endl; 00409 cout << "Track: " << aRecoTrk.id() << " nHit: " 00410 << hitList->nHit() << " nActive: " 00411 << tkFit->nActive() << " chisq/dof: " << 00412 tkFit->chisq()/(tkFit->nActive() - 5) << endl; 00413 TrkRecoTrk& othTrack = trkXRef[trackMost]->track(); 00414 cout << "Track: "<< othTrack.id() << " nHit: " << 00415 othTrack.hits()->nHit() << " nActive: " << 00416 othTrack.fitResult()->nActive() << " chisq/dof: " << 00417 othTrack.fitResult()->chisq() / 00418 (othTrack.fitResult()->nActive() - 5) << endl; 00419 } 00420 00421 if (groupDiff > 0.0) { 00422 // inactivate hits on this track 00423 lRefit = 1; 00424 for (int ii = 0; ii < nMost; ii++) { 00425 TrkHitOnTrk *alink = theseHits[ii]; 00426 TrkFundHit* hit = const_cast<TrkFundHit*> (alink->hit()); 00427 if(hitList->removeHit(hit)){ 00428 alink->hit()->setUnusedHit(alink);//yzhang 2009-11-03 open 00429 } 00430 //alink->setActivity(0); 00431 } 00432 if (8 == tkParam.lPrint) std::cout<<"inactive hits on this track, No."<<aRecoTrk.id()<< std::endl; 00433 } else { 00434 // inactivate hits on other track 00435 refitTrack[trackMost] = 1; 00436 for (int ii = 0; ii < nMost; ii++) { 00437 TrkHitOnTrk *alink = thoseHits[ii]; 00438 TrkFundHit* hit = const_cast<TrkFundHit*> (alink->hit()); 00439 if(hitList->removeHit(hit)){ 00440 alink->hit()->setUnusedHit(alink);//yzhang 2009-11-03 open 00441 } 00442 //alink->setActivity(0); 00443 } 00444 if (8 == tkParam.lPrint) std::cout<<"inactive hits on other track "<< std::endl; 00445 } 00446 delete [] theseHits; 00447 delete [] thoseHits; 00448 00449 } // end if lGroupHits 00450 00451 } // end if nMost > 0 00452 00453 //********* 00454 // Refit this track, if any hits have been dropped 00455 TrkErrCode fitResult; 00456 long index = -1; 00457 idMap.get(aRecoTrk.id(), index); assert (index >= 0); 00458 00459 if (lRefit || refitTrack[index] == 1) { 00460 if (8 == tkParam.lPrint) { 00461 std::cout<<"after group ,refit track"<<aRecoTrk.id()<< std::endl; 00462 } 00463 fitResult = hitList->fit(); 00464 aRecoTrk.status()->addHistory( 00465 TrkErrCode(fitResult.success()?TrkErrCode::succeed:TrkErrCode::fail,14,"Arbitrated"), "MdcTrkRecon"); 00466 if (fitResult.failure() && (8 == tkParam.lPrint )) { 00467 fitResult.print(std::cerr); 00468 } 00469 00470 00471 double chisqperDOF; 00472 bool badFit = true; 00473 if (fitResult.success()) { 00474 badFit = false; 00475 int nDOF = tkFit->nActive() - 5; 00476 if (nDOF > 5){ 00477 chisqperDOF = tkFit->chisq() / nDOF; 00478 }else{ 00479 chisqperDOF = tkFit->chisq(); 00480 } 00481 00482 if (chisqperDOF > tkParam.maxChisq) badFit = true; 00483 if (tkFit->nActive() < tkParam.minHits) badFit = true; 00484 double tem2 = (float) hitList->nHit() - tkFit->nActive(); 00485 if (tkParam.lUseQualCuts) { 00486 if (tem2 >= tkParam.maxNmissTrack) badFit = true; 00487 if (tem2 /float(hitList->nHit()) > tkParam.maxNmissNorm){ 00488 badFit = true; 00489 } 00490 } 00491 if(8== tkParam.lPrint) std::cout<<"fit quality:"<< 00492 " chisqperDof "<<chisqperDOF<<"?>"<<tkParam.maxChisq<< 00493 " nActive "<<tkFit->nActive()<<"?<"<<tkParam.minHits<< 00494 " nHit "<<hitList->nHit()<<" nhit-act "<<tem2<<"?>= nMiss "<<tkParam.maxNmissTrack<< 00495 " hit-act/nhit "<<tem2/float(hitList->nHit())<<"?> MissNorm "<<tkParam.maxNmissNorm 00496 << std::endl; 00497 00498 00499 } 00500 if (8 == tkParam.lPrint) { 00501 cout << "Refitting track " << aRecoTrk.id() << " success = " 00502 << fitResult.success() << "\n"; 00503 } 00504 // If the track no longer passes cuts, delete it 00505 if (fitResult.failure() || badFit ) { 00506 nDeleted++; 00507 // Don't change the track list while we're iterating through it! 00508 // remove(atrack); 00509 //int id = aRecoTrk.id(); 00510 if (8 == tkParam.lPrint) { 00511 cout << "fitResult.failure? "<<fitResult.failure() 00512 <<" badFit? "<<badFit <<" Killing tkNo " << itrack << endl; 00513 } 00514 delete &(atrack->track()); // Delete the RecoTrk inside atrack 00515 atrack->setTrack(0); 00516 trksToKill.push_back(atrack); 00517 continue; 00518 } 00519 } // end if lRefit 00520 00521 if (lGroupHits) goto restart; 00522 00523 } // end loop over tracks 00524 if (8 == tkParam.lPrint) std::cout<<"end of loop over tracks"<< std::endl; 00525 00526 // Remove dead track husks 00527 for (int itk = 0; itk < (int)trksToKill.size(); itk++) { 00528 remove(trksToKill[itk]); 00529 if (8 == tkParam.lPrint) std::cout<<"remode dead track No."<<itk<< std::endl; 00530 } 00531 if (8 == tkParam.lPrint) std::cout<<"---end of arbitrateHits"<< std::endl; 00532 00533 delete [] usedInTrackNum; 00534 delete [] refitTrack; 00535 delete [] trkXRef; 00536 return nDeleted; 00537 }
|
|
Implements MdcTrackListBase. |
|
Implements MdcTrackListBase. 00047 { 00048 // ************************************************************************* 00049 // Forms tracks out of list of superlayer segments 00050 int nTracks = 0; 00051 // Create empty list of stereo segs (to save time) 00052 MdcSegGrouperSt stereoSegs(gm,tkParam.lPrint); 00053 // *** Create r-phi track 00054 00055 // Create list of axial segments, treated as on straight tracks 00056 MdcSegGrouperCsmc straightSegs(gm, tkParam.lPrint); 00057 straightSegs.fillWithSegs(segs); 00058 //std::cout<< "after straight fill " << std::endl; 00059 //segs->plot();//yzhang debug 00060 MdcSeg *seed; 00061 int goodOnly = 1; 00062 MdcTrack* trialTrack = 0; 00063 while (1) { 00064 seed = segs->getSeed(0,goodOnly); 00065 if (seed == 0 && goodOnly == 1) { 00066 segs->resetSeed(gm); 00067 goodOnly = 0; 00068 continue; 00069 } 00070 else if (seed == 0 && goodOnly == 0) { 00071 break; 00072 } 00073 delete trialTrack; 00074 trialTrack = 0; 00075 int success = straightSegs.combineSegs(trialTrack, seed, 00076 context, t0, tkParam.maxSegChisqO); 00077 if (!success) { 00078 //std::cout<< " MdcTrackListCsmc::combineSegs not successed!" << std::endl; 00079 continue; 00080 } 00081 00082 if (tkParam.lPrint > 1) { 00083 trialTrack->track().printAll(cout); 00084 } 00085 // *** End r-phi track-finding 00086 00087 // *** Add stereo to taste 00088 stereoSegs.fillWithSegs(segs, trialTrack); 00089 00090 success = stereoSegs.combineSegs(trialTrack, 0, context, t0, 00091 tkParam.maxSegChisqO); 00092 if (success) { 00093 // Finish 3-d track; 00094 success = finish3d(trialTrack->track()); 00095 // success = finish3d(trialTrack->track()); // GSciolla: try to repeat 00096 } 00097 if (!success) { 00098 //std::cout<< " MdcTrackListCsmc::finish3d not successed!" << std::endl; 00099 continue; 00100 } 00101 00102 if (tkParam.lPrint > 1) { 00103 trialTrack->track().printAll(cout); 00104 } 00105 00106 nTracks++; 00107 append(trialTrack); // Add to list of Tracks 00108 trialTrack = 0; 00109 00110 } // end while(1) 00111 delete trialTrack; 00112 00113 00114 //cout << " Number of Tracks found= " << nTracks ; 00115 //cout << " " << endl; 00116 return nTracks; 00117 }
|
|
Reimplemented in MdcTrackList, and MdcTrackList. |
|
Reimplemented in MdcTrackList, and MdcTrackList. |
|
|
|
00120 { 00121 // ************************************************************************ 00122 int success = 0; 00123 int nParamFit = 0 ; 00124 00125 TrkErrCode fitResult; 00126 00127 // ------------------------ 00128 // 4 param fit (line) 00129 // ------------------------ 00130 nParamFit = 4; 00131 TrkLineMaker makeFit; 00132 makeFit.changeFit(trk); 00133 makeFit.setFlipAndDrop(trk, true, true); 00134 //setCosmic(&trk); // set hot to cosmics OBSOLETE! REMOVE! 00135 fitResult = trk.hits()->fit(); 00136 makeFit.setFlipAndDrop(trk, false, false); 00137 00138 // --------------------------------------------------------------- 00139 // Are there some TDC that can be replaced with later measurements? ( beginning) 00140 // --------------------------------------------------------------- 00141 // Note to Sciolla: multiple hits commented out for now OK! I will put it back later... 00142 /* if( _flags.useMultipleHits ) { 00143 int NHITS = trk.nHit(); 00144 00145 if (fitResult.success()) { 00146 int ifirst = 0; 00147 for (int ihot = 0; ihot < NHITS ; ihot++) { 00148 00149 MdcHitOnTrack* ahot = (MdcHitOnTrack*)trk.hitAList()[ihot]; 00150 double firstTime = ahot->mdcHit()->rawTime(0); 00151 00152 // Get the list of times for this Digi 00153 const MdcDigi* tmpDigi = ahot->mdcHit()->digi() ; 00154 // int nTDChits = tmpDigi->getTdcTimeAList()->length() ; 00155 int nTDChits = tmpDigi->tmdcits() ; 00156 00157 double newTime=0. ; 00158 00159 for( int il = 1; il<nTDChits ; il++){ 00160 00161 // get the time corresponding to this doca ---> distance to time 00162 newTime = tmpDigi->TdcTime(il); 00163 double newDrift = ahot->layer()->timeToDistance(newTime); 00164 00165 double tmp_doca = ahot->dcaToWire(); 00166 double new_diff = fabs(fabs(tmp_doca)-fabs(newDrift)); 00167 double tmp_drift = ahot->drift(); 00168 double tmp_drift2 = 00169 ahot->layer()->timeToDistance(tmpDigi->TdcTime(0)); 00170 double old_diff = fabs(fabs(tmp_doca)-fabs(tmp_drift2)); 00171 00172 if( ((old_diff-new_diff)>0.1) // new time better then old of at least 1 mm 00173 && (new_diff<=0.8) ) { // and newresid < 2 mm 00174 // if(_flags.debug) { 00175 if(1) { 00176 00177 if(ifirst==0) { 00178 ifirst = 1 ; 00179 cout <<" nHits | time1 |fittim | time2 | doca(mm) | drift1 " 00180 <<" drift1(rec) | drift2 | layer | wire | isActive" 00181 << " | chi2 | nactive" << endl; 00182 } 00183 00184 int isAct = 1 ; 00185 if (!ahot->isActive()) { 00186 isAct = 0 ; 00187 // cout<< " >> " ; 00188 // cout << "This hit is not active " << endl ; 00189 } 00190 00191 cout << " " << nTDChits 00192 << " " << tmpDigi->TdcTime(0) 00193 << " " << ahot->fitTime() 00194 << " " << newTime 00195 << " " << tmp_doca*10 00196 << " " << tmp_drift*10 00197 << " " << tmp_drift2*10 00198 << " " << newDrift*10 00199 << " " << ahot->layer()->layNum() 00200 << " " << ahot->wire() 00201 << " " << isAct 00202 << " " << trk.chisq() 00203 << " " << trk.nActive() 00204 << endl; 00205 00206 cout << " old/new diff (mm) = " << old_diff*10 00207 << " | " << new_diff*10 00208 << endl ; 00209 } 00210 00211 ahot->setTimeIndex(il) ; 00212 // cout << "new rawtime =" 00213 // << ahot->mdcHit()->rawTime(0)<< endl ; // check the new value 00214 00215 00216 // store results in ntuple to see improvements ... before ... 00217 HepTuple* tupleMultHits = _manager->ntuple("multiHits"); 00218 00219 tupleMultHits->column("ch2Dof1", trk.chisq()/(trk.nActive() - 4)); 00220 tupleMultHits->column("nActiv1", trk.nActive()); 00221 tupleMultHits->column("resid1",ahot->resid()); 00222 tupleMultHits->column("doca1",ahot->resid()+ahot->drift()); 00223 tupleMultHits->column("time1", ahot->fitTime() ); 00224 tupleMultHits->column("rawtime1", ahot->mdcHit()->rawTime(il) ); 00225 00226 // - make the hit active and usable 00227 ahot->setActivity(true); 00228 00229 // now refit the track ... 00230 fitResult = trk.fit(); 00231 00232 // store results in ntuple to see improvements ... after ... 00233 tupleMultHits->column("ch2Dof2", trk.chisq()/(trk.nActive() - 4)); 00234 tupleMultHits->column("nActiv2", trk.nActive()); 00235 tupleMultHits->column("resid2",ahot->resid()); 00236 tupleMultHits->column("doca2",ahot->resid()+ahot->drift()); 00237 tupleMultHits->column("time2", ahot->fitTime() ); 00238 tupleMultHits->column("rawtime2", ahot->mdcHit()->rawTime(il) ); 00239 tupleMultHits->dumpData(); 00240 } 00241 } 00242 } 00243 } 00244 } // Are there some TDC that can be replaced with later measurements? (end) 00245 */ 00246 00247 const TrkFit* tkFit = trk.fitResult(); 00248 double chisqperDOF = 0.; 00249 if (fitResult.success()) { 00250 int nDOF = tkFit->nActive() - nParamFit; 00251 if (nDOF > nParamFit) 00252 chisqperDOF = tkFit->chisq() / nDOF; 00253 else 00254 chisqperDOF = tkFit->chisq(); 00255 int goodMatch = 1; 00256 if (chisqperDOF > tkParam.maxChisq) goodMatch = 0; 00257 if (tkFit->nActive() < tkParam.minHits) goodMatch = 0; 00258 00259 if (goodMatch) success = 1; 00260 } 00261 00262 return success; 00263 }
|
|
|
|
00541 { 00542 //************************************************************************** 00543 tkParam = tkPar; 00544 }
|
|
00039 {return length();}
|
|
00039 {return length();}
|
|
|
|
00092 { 00093 //************************************************************************* 00094 std::cout<< "nTrack "<<nTrack() << std::endl;//yzhang debug 00095 for (int itrack = 0; itrack < nTrack(); itrack++) { 00096 MdcTrack *atrack = (*this)[itrack]; 00097 if (atrack == NULL) continue; 00098 atrack->track().printAll(cout); 00099 } 00100 }
|
|
Reimplemented from MdcTrackListBase. |
|
Reimplemented from MdcTrackListBase. |
|
00057 {m_d0Cut = d0Cut;}//yzhang add
|
|
00057 {m_d0Cut = d0Cut;}//yzhang add
|
|
|
|
|
|
00059 {m_ptCut = ptCut;}//yzhang add 2009-10-27
|
|
00059 {m_ptCut = ptCut;}//yzhang add 2009-10-27
|
|
00058 {m_z0Cut = z0Cut;}//yzhang add 2010-05-21
|
|
00058 {m_z0Cut = z0Cut;}//yzhang add 2010-05-21
|
|
|
|
00075 { 00076 // *********************************************************************** 00077 int trackId = 0; 00078 for (int itrack = 0; itrack < nTrack(); itrack++) { 00079 MdcTrack* track = (*this)[itrack]; 00080 //tkStat: 0,PatRec 1,MdcxReco 2,Tsf 3,CurlFinder 00081 int tkStat = 0; 00082 track->storeTrack(trackId, trackList, hitList, tkStat); 00083 ++trackId; 00084 } 00085 HepAListDeleteAll(*this); // Discard the husks 00086 removeAll(); 00087 return; 00088 }
|
|
|
|
|
|
|
|
|