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 #include "MdcTrkRecon/MdcMergeDups.h"
00029 #include "GaudiKernel/MsgStream.h"
00030 #include "GaudiKernel/AlgFactory.h"
00031 #include "GaudiKernel/ISvcLocator.h"
00032 #include "GaudiKernel/DataSvc.h"
00033 #include "GaudiKernel/SmartDataPtr.h"
00034 #include "GaudiKernel/IDataProviderSvc.h"
00035 #include "GaudiKernel/IDataManagerSvc.h"
00036 #include "GaudiKernel/PropertyMgr.h"
00037 #include "EventModel/EventHeader.h"
00038 #include "MagneticField/IMagneticFieldSvc.h"
00039 #include "Identifier/MdcID.h"
00040 #include "Identifier/Identifier.h"
00041 #include "MdcRecEvent/RecMdcHit.h"
00042 #include "MdcGeom/Constants.h"
00043 #include "TrkBase/TrkExchangePar.h"
00044 #include "TrkBase/TrkRecoTrk.h"
00045 #include "TrkBase/TrkHitList.h"
00046 #include "TrkBase/TrkErrCode.h"
00047 #include "TrkBase/TrkFit.h"
00048 #include "MdcData/MdcHitMap.h"
00049 #include "MdcData/MdcHitUse.h"
00050 #include "MdcData/MdcRecoHitOnTrack.h"
00051 #include "MdcData/MdcHitOnTrack.h"
00052 #include "TrkFitter/TrkHelixMaker.h"
00053 #include "TrkFitter/TrkContextEv.h"
00054 #include "EvTimeEvent/RecEsTime.h"
00055 #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
00056 #include "MdcTrkRecon/MdcTrack.h"
00057
00058
00059 MdcMergeDups::MdcMergeDups(const std::string& name, ISvcLocator* pSvcLocator) :
00060 Algorithm(name, pSvcLocator)
00061 {
00062 declareProperty("debug", m_debug = 0);
00063
00064 declareProperty("maxDd0InMerge", m_maxDd0InMerge = 2.7);
00065 declareProperty("maxDphi0InMerge", m_maxDphi0InMerge = 0.15);
00066 declareProperty("maxDPdradInMerge", m_maxPdradInMerge= 0.22);
00067 declareProperty("maxRcsInMerge", m_maxRcsInMerge = 18.);
00068
00069 declareProperty("mergePt", m_mergePt = 0.13);
00070 declareProperty("mergeLoadAx", m_mergeLoadAx = 3.);
00071 declareProperty("mergeLoadSt", m_mergeLoadSt = 4.);
00072 declareProperty("mergeOverlapRatio", m_mergeOverlapRatio = 0.7);
00073 }
00074
00075
00076
00077
00078 MdcMergeDups::~MdcMergeDups() {
00079 delete m_bfield;
00080 }
00081
00082 StatusCode MdcMergeDups::beginRun(){
00083
00084 m_gm = MdcDetector::instance();
00085 if(NULL == m_gm) return StatusCode::FAILURE;
00086 return StatusCode::SUCCESS;
00087 }
00088
00089
00090
00091
00092 StatusCode MdcMergeDups::initialize(){
00093 MsgStream log(msgSvc(), name());
00094 log << MSG::INFO << "in initialize()" << endreq;
00095 StatusCode sc;
00096
00097
00098
00099 IMagneticFieldSvc* m_pIMF;
00100 sc = service ("MagneticFieldSvc",m_pIMF);
00101 if(sc != StatusCode::SUCCESS) {
00102 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00103 return StatusCode::FAILURE;
00104 }
00105 m_bfield = new BField(m_pIMF);
00106
00107 return StatusCode::SUCCESS;
00108 }
00109
00110
00111 StatusCode MdcMergeDups::execute() {
00112 MsgStream log(msgSvc(), name());
00113 log << MSG::INFO << "in execute()" << endreq;
00114 setFilterPassed(false);
00115
00116 m_bunchT0 = -999.;
00117 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00118 if (!aevtimeCol || aevtimeCol->size()==0) {
00119 log << MSG::WARNING<< " Could not find RecEsTimeCol"<< endreq;
00120 return StatusCode::SUCCESS;
00121 }
00122
00123 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00124 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00125 m_bunchT0 = (*iter_evt)->getTest();
00126 }
00127
00128
00129 int nMerged = mergeCurl();
00130
00131 if(m_debug>0) {
00132 std::cout<<name()<<": Merged "<<nMerged << " track "<< std::endl;
00133 dumpRecMdcTrack();
00134 }
00135
00136 return StatusCode::SUCCESS;
00137 }
00138
00139
00140 StatusCode MdcMergeDups::finalize(){
00141 MsgStream log(msgSvc(), name());
00142 log << MSG::INFO << "in finalize()" << endreq;
00143
00144 return StatusCode::SUCCESS;
00145 }
00146
00147
00148
00149 int MdcMergeDups::mergeCurl(){
00150
00151
00152 SmartDataPtr<RecMdcTrackCol> trackList(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00153 if (!trackList) return -1;
00154
00155 int needMerge = 0;
00156
00157
00158 RecMdcTrackCol::iterator iterRefTk = trackList->begin();
00159 for (; iterRefTk != trackList->end(); iterRefTk++) {
00160 RecMdcTrack* refTk = *iterRefTk;
00161 if (refTk->stat()<0) continue;
00162 std::vector<RecMdcTrack*> mergeTkList;
00163 mergeTkList.push_back(refTk);
00164
00165
00166 bool curl = false;
00167 int sameParm = 0;
00168 RecMdcTrackCol::iterator iterTestTk = trackList->begin();
00169 for (; iterTestTk != trackList->end(); iterTestTk++) {
00170 RecMdcTrack* testTk = *iterTestTk;
00171 if (iterRefTk == iterTestTk || (testTk->stat()<0)) continue;
00172
00173
00174 if (testByOverlapHit(refTk,testTk)){
00175 if(m_debug>0)std::cout<<__FILE__<<" overlape tk:" <<refTk->trackId()<<" with "<<testTk->trackId()<< std::endl;
00176 mergeTkList.push_back(testTk);
00177 curl = true;
00178 }
00179 sameParm = testByParam(refTk,testTk);
00180 if(sameParm >0) {
00181 if(m_debug>0) std::cout<<__FILE__<<" same param tk:" <<refTk->trackId()<<" with "<<testTk->trackId()<< std::endl;
00182 mergeTkList.push_back(testTk);
00183 }
00184 }
00185 if (mergeTkList.size()>1 && curl) needMerge = doMergeCurl(mergeTkList);
00186 if ((needMerge < 999) && mergeTkList.size()>1 ) needMerge = doMergeLong(mergeTkList);
00187
00188 }
00189
00190
00191 if( needMerge <=0 ) return 0;
00192
00193
00194 iterRefTk = trackList->begin();
00195 int iTk=0;
00196 int nDeleted = 0;
00197 for (; iterRefTk != trackList->end(); ) {
00198 if ( (*iterRefTk)->stat() >= 0 ){
00199 (*iterRefTk)->setTrackId(iTk);
00200 iterRefTk++;
00201 iTk++;
00202 }else {
00203 int id = (*iterRefTk)->trackId();
00204 bool erased = eraseTdsTrack(iterRefTk);
00205 if ( erased ){
00206 nDeleted++;
00207 if(m_debug>0)std::cout<<__FILE__<<" erase track No."<<id<< std::endl;
00208 }else {
00209 if(m_debug>0)std::cout<<__FILE__<<" erase failed !"<< std::endl;
00210 }
00211 }
00212
00213 }
00214 if(m_debug>0) std::cout<<__FILE__<<" After merge save "<<iTk<<" tracks"<< std::endl;
00215
00216 return nDeleted;
00217 }
00218
00219
00220
00221
00222 int MdcMergeDups::testByParam(RecMdcTrack* refTk, RecMdcTrack* testTk){
00223
00224 int overlaped = 0;
00225
00226
00227
00228 double Bz = m_bfield->bFieldZ();
00229 double omega1 = (Constants::c * Bz*refTk->helix(2))/10000.;
00230 double omega2 = (Constants::c * Bz*testTk->helix(2))/10000.;
00231
00232 double phi01 = refTk->helix(1)+Constants::pi/2.;
00233 double phi02 = testTk->helix(1)+Constants::pi/2.;
00234 while(phi01>Constants::twoPi) phi01 -= Constants::twoPi;
00235 while(phi02>Constants::twoPi) phi02 -= Constants::twoPi;
00236 double d01 = -refTk->helix(0);
00237 double d02 = -testTk->helix(0);
00238 double dphi0 = fabs(phi01 - phi02);
00239 double dd0 = fabs(d01 - d02);
00240 double prodo = omega1*omega2;
00241 double r1=100000.;
00242 double r2=100000.;
00243 if (fabs(omega1)>0.00001) r1 = 1.0/fabs(omega1);
00244 if (fabs(omega2)>0.00001) r2 = 1.0/fabs(omega2);
00245 double pdrad = fabs((r1-r2)/(r1+r2)) ;
00246
00247 if (2==m_debug){
00248 std::cout << " fabs(d01 - d02) " << fabs(d01 - d02) << std::endl;
00249 std::cout << " fabs(phi01-phi02) " << fabs(phi01-phi02) << std::endl;
00250 }
00251
00252 if ( (prodo > 0.) && (dd0 < m_maxDd0InMerge) && (dphi0 < m_maxDphi0InMerge) &&
00253 (pdrad < m_maxPdradInMerge)) {
00254 overlaped = 1;
00255 }
00256
00257
00258 if ( (prodo < 0.) && (fabs(d01+d02) < 4.0) && (dd0 > 47.0) &&
00259 (fabs( dphi0 - Constants::pi) < m_maxDphi0InMerge)
00260 && (pdrad < m_maxPdradInMerge)) {
00261 overlaped = 2;
00262 }
00263
00264 return overlaped;
00265 }
00266
00267
00268
00269
00270 int MdcMergeDups::doMergeLong(std::vector<RecMdcTrack*> mergeTkList){
00271
00272
00273 double minRcs=999.;
00274 int bestTkId=999;
00275 RecMdcTrack* bestTk=NULL;
00276 std::vector<RecMdcTrack*>::iterator itTk = mergeTkList.begin();
00277 for (int iTk=0; itTk != mergeTkList.end(); itTk++,iTk++){
00278 RecMdcTrack* tk = (*itTk);
00279 double chi2 = tk->chi2();
00280 double ndf = tk->ndof();
00281 if(chi2/ndf < minRcs) {
00282 bestTkId = tk->trackId();
00283 bestTk = tk;
00284 }
00285 }
00286 if (minRcs < m_maxRcsInMerge) return bestTkId;
00287 bestTk->setStat(-1);
00288
00289 return 999;
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 int MdcMergeDups::testByOverlapHit(RecMdcTrack* refTk, RecMdcTrack* testTk){
00349
00350 int overlaped = 0;
00351 if ((testTk->pxy() >= m_mergePt) || (refTk->pxy() >= m_mergePt)) return overlaped;
00352
00353 HitRefVec testHits = testTk->getVecHits();
00354 int nHit = testHits.size();
00355 int nOverlap = 0;
00356
00357 HitRefVec::iterator iterHit = testHits.begin();
00358 for (; iterHit != testHits.end(); iterHit++) {
00359 RecMdcHit* hit = *iterHit;
00360
00361
00362 double load = m_mergeLoadAx;
00363 bool isStLayer = (m_gm->Layer(MdcID::layer(hit->getMdcId()))->view() != 0);
00364 if(isStLayer) load = m_mergeLoadSt;
00365
00366
00367 double vx0 = refTk->getVX0();
00368 double vy0 = refTk->getVY0();
00369 double vz0 = refTk->getVZ0();
00370 double dr = refTk->helix(0);
00371 double phi0 = refTk->helix(1);
00372 double Bz = m_bfield->bFieldZ();
00373 double r = 10000./ (Constants::c * Bz*refTk->helix(2));
00374 double dz = refTk->helix(3);
00375 double tanl = refTk->helix(4);
00376
00377
00378 double xc = vx0 + (dr + r) * cos(phi0);
00379 double yc = vy0 + (dr + r) * sin(phi0);
00380
00381
00382 double zHit = hit->getZhit();
00383 double phi = (vz0 + dz - zHit) / (r * tanl);
00384 double xHit = vx0 + dr*cos(phi0) + r*(cos(phi0) - cos(phi0+phi));
00385 double yHit = vy0 + dr*sin(phi0) + r*(sin(phi0) - sin(phi0+phi));
00386
00387
00388 double dx = xc - xHit;
00389 double dy = yc - yHit;
00390 double dHit2Center = sqrt(dx * dx + dy * dy);
00391 double rTk = fabs(r);
00392
00393
00394 if ( (dHit2Center>(rTk - load)) && (dHit2Center<(rTk + load))) nOverlap++;
00395 }
00396
00397 if ( nOverlap<=0 ) return overlaped;
00398
00399 double overlapRatio = double(nOverlap) / double(nHit);
00400
00401 if (overlapRatio > m_mergeOverlapRatio) overlaped = 1;
00402
00403 return overlaped;
00404 }
00405
00406
00407 int MdcMergeDups::doMergeCurl(std::vector<RecMdcTrack*> mergeTkList){
00408
00409 int innerMostTkId = 999;
00410 RecMdcTrack* innerMostTk = NULL;
00411 unsigned innerMostLayerOfTk = 999;
00412 std::vector<RecMdcTrack*>::iterator itTk = mergeTkList.begin();
00413 for (int iTk=0; itTk != mergeTkList.end(); itTk++,iTk++){
00414 RecMdcTrack* tk = (*itTk);
00415 unsigned innerMostLayer = 999;
00416 for (unsigned iHit = 0; iHit < tk->getVecHits().size(); iHit++) {
00417 unsigned layer = MdcID::layer(tk->getVecHits()[iHit]->getMdcId());
00418 if (layer < innerMostLayer) innerMostLayer=layer;
00419 }
00420
00421 if(m_debug>0)std::cout<<__FILE__<<" to be merged track id="<<tk->trackId()<< std::endl;
00422
00423 if(innerMostLayer < innerMostLayerOfTk){
00424 innerMostTkId = iTk;
00425 innerMostTk = tk;
00426 }else if (innerMostLayer == innerMostLayerOfTk) {
00427
00428 if (tk->helix(3) < innerMostTk->helix(3)){
00429 innerMostTkId = iTk;
00430 innerMostTk = tk;
00431 }
00432 }
00433 }
00434 innerMostTk->setStat(-1);
00435
00436 return innerMostTkId;
00437 }
00438
00439 void MdcMergeDups::store(TrkRecoTrk* aTrack){
00440 SmartDataPtr<RecMdcTrackCol> trackList(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00441 if (!trackList) return;
00442 SmartDataPtr<RecMdcHitCol> hitList(eventSvc(),EventModel::Recon::RecMdcHitCol);
00443 if (!hitList) return;
00444
00445 assert (aTrack != NULL);
00446 TrkExchangePar helix = aTrack->fitResult()->helix(0.);
00447
00448 if(m_debug>1)std::cout<<__FILE__<<" STORED"<< std::endl;
00449 MdcTrack mdcTrack(aTrack);
00450
00451 int tkStat = 4;
00452 mdcTrack.storeTrack(-1, trackList, hitList, tkStat);
00453 }
00454
00455 bool MdcMergeDups::eraseTdsTrack(RecMdcTrackCol::iterator tk){
00456 SmartDataPtr<RecMdcTrackCol> trackList(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00457 if (!trackList) return false;
00458 SmartDataPtr<RecMdcHitCol> hitList(eventSvc(),EventModel::Recon::RecMdcHitCol);
00459 if (!hitList) return false;
00460 HitRefVec hits = (*tk)->getVecHits();
00461 HitRefVec::iterator iterHit = hits.begin();
00462 for (; iterHit != hits.end(); iterHit++) {
00463
00464 }
00465 trackList->erase(tk);
00466 return true;
00467 }
00468
00469
00470 void MdcMergeDups::dumpRecMdcTrack(){
00471 SmartDataPtr<RecMdcTrackCol> trackList(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00472 if (!trackList) return;
00473 if (trackList->size() != 4 ) setFilterPassed(true);
00474 std::cout<<"N track after Merged = "<<trackList->size() << std::endl;
00475 if (m_debug <=1) return;
00476 RecMdcTrackCol::iterator it = trackList->begin();
00477 for (;it!= trackList->end();it++){
00478 RecMdcTrack *tk = *it;
00479 std::cout<< "//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
00480 cout <<" d0 "<<tk->helix(0)
00481 <<" phi0 "<<tk->helix(1)
00482 <<" cpa "<<tk->helix(2)
00483 <<" z0 "<<tk->helix(3)
00484 <<" tanl "<<tk->helix(4)
00485 <<endl;
00486 std::cout<<" q "<<tk->charge()
00487 <<" theta "<<tk->theta()
00488 <<" phi "<<tk->phi()
00489 <<" x0 "<<tk->x()
00490 <<" y0 "<<tk->y()
00491 <<" z0 "<<tk->z()
00492 <<" r0 "<<tk->r()
00493 <<endl;
00494 std::cout <<" p "<<tk->p()
00495 <<" pt "<<tk->pxy()
00496 <<" px "<<tk->px()
00497 <<" py "<<tk->py()
00498 <<" pz "<<tk->pz()
00499 <<endl;
00500 std::cout<<" tkStat "<<tk->stat()
00501 <<" chi2 "<<tk->chi2()
00502 <<" ndof "<<tk->ndof()
00503 <<" nhit "<<tk->getNhits()
00504 <<" nst "<<tk->nster()
00505 <<endl;
00506
00507
00508
00509
00510 int nhits = tk->getVecHits().size();
00511 std::cout<<nhits <<" Hits: " << std::endl;
00512 for(int ii=0; ii <nhits ; ii++){
00513 Identifier id(tk->getVecHits()[ii]->getMdcId());
00514 int layer = MdcID::layer(id);
00515 int wire = MdcID::wire(id);
00516 cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
00517 <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
00518 }
00519 std::cout << " "<< std::endl;
00520 }
00521 std::cout << " "<< std::endl;
00522 }