00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <math.h>
00014 #include <iostream>
00015 #include <vector>
00016 #include "MdcGeom/Constants.h"
00017 #include "MdcGeom/BesAngle.h"
00018 #include "MdcData/MdcHit.h"
00019 #include "MdcData/MdcHitOnTrack.h"
00020 #include "MdcGeom/MdcDetector.h"
00021 #include "MdcTrkRecon/MdcTrackListBase.h"
00022 #include "MdcTrkRecon/MdcMap.h"
00023 #include "MdcTrkRecon/MdcTrack.h"
00024 #include "TrkBase/TrkErrCode.h"
00025 #include "TrkBase/TrkRecoTrk.h"
00026 #include "TrkBase/TrkFit.h"
00027 #include "TrkBase/TrkFitStatus.h"
00028 #include "TrkBase/TrkHitList.h"
00029 #include "TrkBase/TrkExchangePar.h"
00030 #include "MdcRecoUtil/Pdt.h"
00031 #include "MdcGeom/BesAngle.h"
00032 #include "MdcData/MdcRecoHitOnTrack.h"
00033 #include "MdcRawEvent/MdcDigi.h"
00034 #include "CLHEP/Matrix/SymMatrix.h"
00035 #include "CLHEP/Vector/ThreeVector.h"
00036 #include "TrkFitter/TrkContextEv.h"
00037 #include "TrkBase/TrkRep.h"
00038 #include "Identifier/MdcID.h"
00039 #include "Identifier/Identifier.h"
00040 #include "AIDA/IHistogram1D.h"
00041 #include "AIDA/IHistogram2D.h"
00042 #include "CLHEP/Geometry/Point3D.h"
00043 #include "CLHEP/Geometry/Vector3D.h"
00044 #include "TrkBase/TrkDifTraj.h"
00045
00046 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00047
00048 typedef HepGeom::Point3D<double> HepPoint3D;
00049 #endif
00050
00051 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00052
00053 typedef HepGeom::Vector3D<double> HepVector3D;
00054 #endif
00055
00056 double MdcTrackListBase::m_d0Cut = -999.;
00057 double MdcTrackListBase::m_z0Cut = -999.;
00058 double MdcTrackListBase::m_ptCut = -999.;
00059
00060 #include "GaudiKernel/NTuple.h"
00061
00062
00063 MdcTrackListBase::MdcTrackListBase(const MdcTrackParams &tkPar) {
00064
00065 tkParam = tkPar;
00066 return;
00067 }
00068
00069
00070 MdcTrackListBase::~MdcTrackListBase() {}
00071
00072
00073
00074 void
00075 MdcTrackListBase::store( RecMdcTrackCol *trackList, RecMdcHitCol *hitList){
00076
00077 int trackId = 0;
00078 for (int itrack = 0; itrack < nTrack(); itrack++) {
00079 MdcTrack* track = (*this)[itrack];
00080
00081 int tkStat = 0;
00082 track->storeTrack(trackId, trackList, hitList, tkStat);
00083 ++trackId;
00084 }
00085 HepAListDeleteAll(*this);
00086 removeAll();
00087 return;
00088 }
00089
00090
00091 void
00092 MdcTrackListBase::plot() const {
00093
00094 std::cout<< "nTrack "<<nTrack() << std::endl;
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 }
00101
00102
00103 int
00104 MdcTrackListBase::arbitrateHits() {
00105
00106
00107
00108
00109
00110
00111
00112
00113
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
00126 int* usedInTrackNum = new int [nTrack()];
00127
00128 MdcTrack** trkXRef = new MdcTrack* [nTrack()];
00129
00130 int *refitTrack = new int [nTrack()];
00131 for (int i = 0; i < nTrack(); i++) {
00132 refitTrack[i] = 0;
00133 }
00134
00135
00136 int itrack;
00137 for (itrack = 0; itrack < nTrack(); itrack++) {
00138 MdcTrack *atrack = (*this)[itrack];
00139 if (atrack == 0) continue;
00140 idMap.put(atrack->track().id(), itrack);
00141 trkXRef[itrack] = atrack;
00142 }
00143
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
00160 int nPrev = 0;
00161 int nHitDeleted = 0;
00162 int maxGapLength = 0;
00163 int nGapGE2= 0;
00164 int nGapGE3= 0;
00165 int nHitInLayer[43];
00166 int nDeleteInLayer[43];
00167 for(int i=0;i<43;i++){
00168 nHitInLayer[i]=0;
00169 nDeleteInLayer[i]=0;
00170 }
00171 if(8 == tkParam.lPrint) std::cout<< "--arbitrate--"<<std::endl;
00172 for (TrkHitList::hot_iterator ihit(hitList->begin()); ihit != hitList->end(); ++ihit){
00173 int nUsed = ihit->hit()->nUsedHits();
00174 if (8 == tkParam.lPrint){
00175 std::cout<<"nUsed="<<nUsed<<":";
00176 ihit->hit()->printAll(std::cout);
00177 }
00178 if (8 == tkParam.lPrint) {
00179 double deltaChi = -999;
00180 ihit->getFitStuff(deltaChi);
00181 std::cout<< "deltaChi="<<deltaChi<<std::endl;
00182 }
00183 int layer = ihit->layerNumber();
00184 nHitInLayer[layer]++;
00185
00186 if (!ihit->isActive()) {
00187
00188
00189
00190 if(tkParam.lRemoveInActive ) {
00191 nDeleteInLayer[layer]++;
00192 if (8 == tkParam.lPrint) {
00193 std::cout<< "=remove above inactive "<<std::endl;
00194 }
00195 TrkFundHit* hit = const_cast<TrkFundHit*> (ihit->hit());
00196 hitList->removeHit(hit);
00197 if(ihit == hitList->end()) break;
00198 --ihit;
00199 }
00200 continue;
00201 }
00202 if (nUsed > 1) {
00203 bool wasUsed = false;
00204 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator> q =
00205 ihit->hit()->getUsedHits();
00206 for (TrkFundHit::hot_iterator i = q.first; i != q.second; ++i) {
00207 if ( !i->isActive() ) continue;
00208 TrkRecoTrk * recoTrk=i->parentTrack();
00209 int id = recoTrk->id();
00210 if (id == aRecoTrk.id()) continue;
00211 long index = 0;
00212 idMap.get(id, index);
00213 assert(index >= 0);
00214 usedInTrackNum[index]++;
00215 if (8 == tkParam.lPrint){
00216 std::cout<<" track "<<itrack<<"&" <<index
00217 << " shared hits "<<usedInTrackNum[index]<<":";
00218 ihit->printAll(std::cout);
00219 }
00220 wasUsed = true;
00221 }
00222 if (wasUsed) nPrev++;
00223 }
00224 }
00225
00226 int testGap = 0;
00227
00228 for (int i=0;i<43;i++){
00229
00230 if (8 == tkParam.lPrint) {
00231 std::cout<<i<<" nHitInLayer "<<nHitInLayer[i]
00232 <<" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl;
00233 }
00234
00235 if(nHitInLayer[i]>0 && (nHitInLayer[i]-nDeleteInLayer[i])==0) {
00236
00237 nHitDeleted++;
00238 if (8 == tkParam.lPrint) {
00239 cout << "rec hits have been deleted in this layer"<<std::endl;
00240 }
00241 testGap++;
00242
00243 }else if(nHitInLayer[i]==0){
00244
00245 testGap++;
00246
00247 }else{
00248
00249
00250 if(testGap>=2){
00251 nGapGE2++;
00252 if(testGap>=3){ nGapGE3++; }
00253 if(testGap>maxGapLength) maxGapLength=testGap;
00254
00255 }
00256 testGap=0;
00257 }
00258 }
00259
00260 bool toBeDeleted = false;
00261
00262 if(tkParam.lPrint>1) std::cout<< "arbitrateHits tkNo:"<<itrack<<" nGapGE2= "<<nGapGE2 << " nGapGE3= "<<nGapGE3 << " maxGapLength= "<<maxGapLength<<std::endl;
00263
00264
00265 if (nHitDeleted >= tkParam.nHitDeleted) {
00266 if (tkParam.lPrint>1) {
00267 cout << "arbitrateHits: nHitDeleted "<<nHitDeleted<<" >= "<<tkParam.nHitDeleted
00268 <<" Killing tkNo " << itrack << endl;
00269 }
00270 toBeDeleted = true;
00271 }
00272
00273
00274
00275 if (nGapGE2 >= tkParam.nGapGE2) {
00276 if (tkParam.lPrint>1) {
00277 cout << "arbitrateHits: nGapGE2 "<<nGapGE2<<" >= "<<tkParam.nGapGE2 <<" Killing tkNo " << itrack << endl;
00278 }
00279 toBeDeleted = true;
00280 }
00281 if (nGapGE3 >= tkParam.nGapGE3) {
00282 if (tkParam.lPrint>1) {
00283 cout << "arbitrateHits: nGapGE3 "<<nGapGE3<<" >= "<<tkParam.nGapGE3 <<" Killing tkNo " << itrack << endl;
00284 }
00285 toBeDeleted = true;
00286 }
00287 if (maxGapLength >= tkParam.maxGapLength) {
00288 if (tkParam.lPrint>1) {
00289 cout << "arbitrateHits: maxGapLength "<<maxGapLength<<" >= "<<tkParam.maxGapLength<<" Killing tkNo " << itrack << endl;
00290 }
00291 toBeDeleted = true;
00292 }
00293
00294 if(toBeDeleted){
00295 nDeleted++;
00296 delete &(atrack->track());
00297 atrack->setTrack(0);
00298 trksToKill.push_back(atrack);
00299 continue;
00300 }
00301
00302
00303
00304 int nMost = 0;
00305 int trackMost = 0;
00306 for (int ii = 0; ii < nTrack(); ii++) {
00307 if (8 == tkParam.lPrint){
00308 std::cout<<"tk:"<<itrack<<"&"<<ii
00309 <<" shared "<<usedInTrackNum[ii]<<" hits "<< std::endl;
00310 }
00311 if (usedInTrackNum[ii] > nMost) {
00312 nMost = usedInTrackNum[ii];
00313 trackMost = ii;
00314 }
00315 }
00316
00317
00318 if (trackMost == trackOld) {
00319 std::cout << "ErrMsg(error) MdcTrackListBase:"
00320 << "Something ghastly happened in MdcTrackListBase::arbitrateHits"
00321 << std::endl;
00322 return 0;
00323 }
00324 trackOld = trackMost;
00325
00326
00327
00328
00329 double groupDiff = 0.0;
00330
00331 int nFound = 0;
00332 TrkHitOnTrk **theseHits = 0;
00333 TrkHitOnTrk **thoseHits = 0;
00334 int lGroupHits = 0;
00335
00336 if (nMost >= tkParam.nOverlap) {
00337 if (8 == tkParam.lPrint){
00338 std::cout<<"track "<<trackMost<<" shared "<<nMost<<" hits > Cut nOverlap "
00339 <<tkParam.nOverlap<<", group hits!"<<std::endl;
00340 }
00341 lGroupHits = 1;
00342 theseHits = new TrkHitOnTrk*[nMost];
00343 thoseHits = new TrkHitOnTrk*[nMost];
00344 }
00345
00346
00347
00348
00349
00350 if(8 == tkParam.lPrint) std::cout<<"Go back through hits, looking up overlap hits"<< std::endl;
00351 if (nMost > 0) {
00352 if (8 == tkParam.lPrint) std::cout<<" nHits= "<< hitList->nHit()<< std::endl;
00353 for (TrkHitList::hot_iterator ihit(hitList->begin()); ihit != hitList->end(); ++ihit) {
00354 int nUsed = ihit->hit()->nUsedHits();
00355
00356 if (8 == tkParam.lPrint){
00357 std::cout<< "--hit go back, nUsed="<<nUsed<<":";
00358 ihit->hit()->printAll(std::cout);
00359 }
00360
00361
00362 if (nUsed < 2) { continue; }
00363
00364
00365 if (!ihit->isActive()) {
00366 if (8 == tkParam.lPrint){ std::cout<<"act=0 continue"<<std::endl; }
00367 continue;
00368 }
00369
00370
00371 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator> q = ihit->hit()->getUsedHits();
00372 while (q.first!=q.second) {
00373 int dropThisHit = 0;
00374 TrkHitOnTrk *otherHot = const_cast<TrkHitOnTrk*>((--q.second).get());
00375 TrkRecoTrk *otherTrack = otherHot->parentTrack();
00376
00377 if (!otherHot->isActive()) continue;
00378
00379
00380 if ( &aRecoTrk == otherTrack) continue;
00381 int otherId = otherTrack->id();
00382 long otherIndex = -1;
00383 idMap.get(otherId, otherIndex); assert(otherIndex >= 0);
00384
00385
00386 if (lGroupHits && otherIndex != trackMost) continue;
00387
00388 if (lGroupHits) {
00389 if (8 == tkParam.lPrint) {
00390 std::cout<<"group hits "<< std::endl;
00391 }
00392
00393
00394
00395
00396 int aDof = tkFit->nActive() - 5;
00397 assert (otherTrack->fitResult() != 0);
00398 int otherDof = otherTrack->fitResult()->nActive() - 5;
00399 if (aDof <= 0) {groupDiff = 999;}
00400 else if (otherDof <= 0) {groupDiff = -999;}
00401 else {
00402 groupDiff += ihit->resid(0) * ihit->resid(0) * ihit->weight() /
00403 aDof -
00404 otherHot->resid(0) * otherHot->resid(0) * otherHot->weight() /
00405 otherDof;
00406 }
00407 theseHits[nFound] = const_cast<TrkHitOnTrk*>(ihit.get());
00408 thoseHits[nFound] = otherHot;
00409 nFound++;
00410 dropThisHit = 1;
00411 } else {
00412
00413 if (8 == tkParam.lPrint) {
00414 std::cout<<"handle hits individually"<< std::endl;
00415 }
00416 nFound++;
00417 if (fabs(ihit->resid(0)) > fabs(otherHot->resid(0)) ) {
00418
00419 lRefit = 1;
00420
00421
00422 const_cast<TrkHitOnTrk*>(ihit.get())->setActivity(0);
00423 dropThisHit = 1;
00424 if (8 == tkParam.lPrint) {
00425 std::cout<<"dorp hit ";
00426 const_cast<TrkHitOnTrk*>(ihit.get())->print(std::cout);
00427 }
00428 break;
00429 } else {
00430
00431 refitTrack[otherIndex] = 1;
00432
00433 otherHot->setActivity(0);
00434 if (8 == tkParam.lPrint) {
00435 std::cout<<"inactive hit on other track";
00436 const_cast<TrkHitOnTrk*>(ihit.get())->print(std::cout);
00437 }
00438 break;
00439 }
00440 }
00441
00442 if (dropThisHit == 1) break;
00443
00444 }
00445
00446
00447 if (lGroupHits && nFound == nMost || nFound == nPrev) {
00448 if (8 == tkParam.lPrint) {
00449 std::cout<<"we've found all of the shared hits on this track,Quit"<<std::endl;
00450 }
00451 break;
00452 }
00453
00454 }
00455
00456
00457 if (lGroupHits) {
00458 if (8 == tkParam.lPrint) {
00459 cout << "nGroup: " << nMost << " groupDiff: " << groupDiff << endl;
00460 cout << "Track: " << aRecoTrk.id() << " nHit: "
00461 << hitList->nHit() << " nActive: "
00462 << tkFit->nActive() << " chisq/dof: " <<
00463 tkFit->chisq()/(tkFit->nActive() - 5) << endl;
00464 TrkRecoTrk& othTrack = trkXRef[trackMost]->track();
00465 cout << "Track: "<< othTrack.id() << " nHit: " <<
00466 othTrack.hits()->nHit() << " nActive: " <<
00467 othTrack.fitResult()->nActive() << " chisq/dof: " <<
00468 othTrack.fitResult()->chisq() /
00469 (othTrack.fitResult()->nActive() - 5) << endl;
00470 }
00471
00472 if (groupDiff > 0.0) {
00473
00474 lRefit = 1;
00475 for (int ii = 0; ii < nMost; ii++) {
00476 TrkHitOnTrk *alink = theseHits[ii];
00477 TrkFundHit* hit = const_cast<TrkFundHit*> (alink->hit());
00478 hitList->removeHit(hit);
00479
00480 }
00481 if (8 == tkParam.lPrint) std::cout<<"inactive hits on this track, No."<<aRecoTrk.id()<< std::endl;
00482 } else {
00483
00484 refitTrack[trackMost] = 1;
00485 for (int ii = 0; ii < nMost; ii++) {
00486 TrkHitOnTrk *alink = thoseHits[ii];
00487 TrkFundHit* hit = const_cast<TrkFundHit*> (alink->hit());
00488 hitList->removeHit(hit);
00489
00490 }
00491 if (8 == tkParam.lPrint) std::cout<<"inactive hits on other track "<< std::endl;
00492 }
00493 delete [] theseHits;
00494 delete [] thoseHits;
00495
00496 }
00497
00498 }
00499
00500
00501
00502 TrkErrCode fitResult;
00503 long index = -1;
00504 idMap.get(aRecoTrk.id(), index); assert (index >= 0);
00505
00506 if (lRefit || refitTrack[index] == 1) {
00507 if (8 == tkParam.lPrint) {
00508 std::cout<<"after group ,refit track"<<aRecoTrk.id()<< std::endl;
00509 }
00510 fitResult = hitList->fit();
00511 aRecoTrk.status()->addHistory(
00512 TrkErrCode(fitResult.success()?TrkErrCode::succeed:TrkErrCode::fail,14,"Arbitrated"), "MdcTrkRecon");
00513 if (fitResult.failure() && (8 == tkParam.lPrint )) {
00514 fitResult.print(std::cerr);
00515 }
00516
00517
00518 double chisqperDOF;
00519 bool badFit = true;
00520 if (fitResult.success()) {
00521 badFit = false;
00522 int nDOF = tkFit->nActive() - 5;
00523 if (nDOF > 5){
00524 chisqperDOF = tkFit->chisq() / nDOF;
00525 }else{
00526 chisqperDOF = tkFit->chisq();
00527 }
00528
00529 if (chisqperDOF > tkParam.maxChisq) badFit = true;
00530 if (tkFit->nActive() < tkParam.minHits) badFit = true;
00531 double tem2 = (float) hitList->nHit() - tkFit->nActive();
00532 if (tkParam.lUseQualCuts) {
00533 if (tem2 >= tkParam.maxNmissTrack) badFit = true;
00534 if (tem2 /float(hitList->nHit()) > tkParam.maxNmissNorm){
00535 badFit = true;
00536 }
00537 }
00538 if(8== tkParam.lPrint) std::cout<<"fit quality:"<<
00539 " chisqperDof "<<chisqperDOF<<"?>"<<tkParam.maxChisq<<
00540 " nActive "<<tkFit->nActive()<<"?<"<<tkParam.minHits<<
00541 " nHit "<<hitList->nHit()<<" nhit-act "<<tem2<<"?>= nMiss "<<tkParam.maxNmissTrack<<
00542 " hit-act/nhit "<<tem2/float(hitList->nHit())<<"?> MissNorm "<<tkParam.maxNmissNorm
00543 << std::endl;
00544
00545
00546 }
00547 if (8 == tkParam.lPrint) {
00548 cout << "Refitting track " << aRecoTrk.id() << " success = "
00549 << fitResult.success() << "\n";
00550 }
00551
00552 if (fitResult.failure() || badFit ) {
00553 nDeleted++;
00554
00555
00556
00557 if (8 == tkParam.lPrint) {
00558 cout << "fitResult.failure? "<<fitResult.failure()
00559 <<" badFit? "<<badFit <<" Killing tkNo " << itrack << endl;
00560 }
00561 delete &(atrack->track());
00562 atrack->setTrack(0);
00563 trksToKill.push_back(atrack);
00564 continue;
00565 }
00566 }
00567
00568 if (lGroupHits) goto restart;
00569
00570 }
00571 if (8 == tkParam.lPrint) std::cout<<"end of loop over tracks"<< std::endl;
00572
00573
00574 for (int itk = 0; itk < (int)trksToKill.size(); itk++) {
00575 remove(trksToKill[itk]);
00576 if (8 == tkParam.lPrint) std::cout<<"remode dead track No."<<itk<< std::endl;
00577 }
00578 if (8 == tkParam.lPrint) std::cout<<"---end of arbitrateHits"<< std::endl;
00579
00580 delete [] usedInTrackNum;
00581 delete [] refitTrack;
00582 delete [] trkXRef;
00583 return nDeleted;
00584 }
00585
00586
00587 void
00588 MdcTrackListBase::newParams(const MdcTrackParams &tkPar) {
00589
00590 tkParam = tkPar;
00591 }
00592
00593
00594 void
00595 MdcTrackListBase::remove( MdcTrack *atrack ) {
00596
00597 if (atrack != 0) {
00598 HepAList<MdcTrack>::remove( atrack );
00599 delete atrack;
00600 }
00601 }
00602
00603
00604 void
00605 MdcTrackListBase::transferTrack() {
00606
00607
00608 }