Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

MdcMergeDups Class Reference

#include <MdcMergeDups.h>

List of all members.

Public Member Functions

StatusCode beginRun ()
StatusCode beginRun ()
int doMergeCurl (std::vector< RecMdcTrack * > mergeTkList)
int doMergeCurl (std::vector< RecMdcTrack * > mergeTkList)
int doMergeLong (std::vector< RecMdcTrack * > mergeTkList)
int doMergeLong (std::vector< RecMdcTrack * > mergeTkList)
void dumpRecMdcTrack ()
void dumpRecMdcTrack ()
bool eraseTdsTrack (RecMdcTrackCol::iterator tk)
bool eraseTdsTrack (RecMdcTrackCol::iterator tk)
StatusCode execute ()
StatusCode execute ()
StatusCode finalize ()
StatusCode finalize ()
StatusCode initialize ()
StatusCode initialize ()
 MdcMergeDups (const std::string &name, ISvcLocator *pSvcLocator)
 MdcMergeDups (const std::string &name, ISvcLocator *pSvcLocator)
int mergeCurl (void)
int mergeCurl (void)
int mergeDups (void)
int mergeDups (void)
void store (TrkRecoTrk *aTrack)
void store (TrkRecoTrk *aTrack)
int testByOverlapHit (RecMdcTrack *refTk, RecMdcTrack *testTk)
int testByOverlapHit (RecMdcTrack *refTk, RecMdcTrack *testTk)
int testByParam (RecMdcTrack *refTk, RecMdcTrack *testTk)
int testByParam (RecMdcTrack *refTk, RecMdcTrack *testTk)
virtual ~MdcMergeDups ()
virtual ~MdcMergeDups ()

Private Attributes

BField * m_bfield
BField * m_bfield
double m_bunchT0
double m_Bz
int m_debug
const MdcDetector * m_gm
const MdcDetector * m_gm
double m_maxDd0InMerge
double m_maxDphi0InMerge
double m_maxPdradInMerge
double m_maxRcsInMerge
double m_mergeLoadAx
double m_mergeLoadSt
double m_mergeOverlapRatio
double m_mergePt

Constructor & Destructor Documentation

MdcMergeDups::MdcMergeDups (  const std::string &  name,
ISvcLocator *  pSvcLocator
) 
  

00059                                                                           :
00060   Algorithm(name, pSvcLocator)
00061 {
00062   declareProperty("debug",              m_debug = 0);
00063   //cuts for mergeDups()
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   //cuts for mergeCurl()
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 }

MdcMergeDups::~MdcMergeDups (   )  [virtual]
  

00078                             {
00079   delete m_bfield;
00080 }

MdcMergeDups::MdcMergeDups (  const std::string &  name,
ISvcLocator *  pSvcLocator
) 
 

virtual MdcMergeDups::~MdcMergeDups (   )  [virtual]
 

Member Function Documentation

StatusCode MdcMergeDups::beginRun (   ) 
 

StatusCode MdcMergeDups::beginRun (   ) 
  

00082                                  {  
00083   //Detector geometry 
00084   m_gm = MdcDetector::instance();
00085   if(NULL == m_gm) return StatusCode::FAILURE;
00086   return StatusCode::SUCCESS;
00087 }

int MdcMergeDups::doMergeCurl (  std::vector< RecMdcTrack * >  mergeTkList  ) 
 

int MdcMergeDups::doMergeCurl (  std::vector< RecMdcTrack * >  mergeTkList  ) 
  

00407                                                                 {
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     // test inner most layer id; if same, test dz
00423     if(innerMostLayer < innerMostLayerOfTk){
00424       innerMostTkId = iTk;
00425       innerMostTk = tk;
00426     }else if (innerMostLayer == innerMostLayerOfTk) {
00427       // test by dz
00428       if (tk->helix(3) < innerMostTk->helix(3)){
00429         innerMostTkId = iTk;
00430         innerMostTk = tk;
00431       }
00432     }
00433   }//end of for mergeTkList
00434   innerMostTk->setStat(-1);
00435 
00436   return innerMostTkId;
00437 }

int MdcMergeDups::doMergeLong (  std::vector< RecMdcTrack * >  mergeTkList  ) 
 

int MdcMergeDups::doMergeLong (  std::vector< RecMdcTrack * >  mergeTkList  ) 
  

00271                                                                 {
00272   //-------------------------------------------------------------------
00273   //merge hitlist
00274   double minRcs=999.;
00275   int bestTkId=999;
00276   RecMdcTrack* bestTk=NULL;
00277   std::vector<RecMdcTrack*>::iterator itTk = mergeTkList.begin();
00278   for (int iTk=0; itTk != mergeTkList.end(); itTk++,iTk++){
00279     RecMdcTrack* tk = (*itTk);
00280     double chi2 = tk->chi2();
00281     double ndf  = tk->ndof();
00282     if(chi2/ndf < minRcs) {
00283       bestTkId = tk->trackId(); 
00284       bestTk = tk;
00285     }
00286   }
00287   if (minRcs < m_maxRcsInMerge) return bestTkId;
00288   bestTk->setStat(-1);
00289 
00290   return 999;
00291   //FIXME
00292   /*
00293   //fit with track parameter respectively
00294   MdcxFittedHel fit1(dcxhlist, *iptr);
00295   MdcxFittedHel fit2(dcxhlist, *trkl[j]); 
00296   int uf = 0; 
00297   //get a best fit 
00298   if ( !fit1.Fail() && (fit1.Rcs()<m_maxRcsInMerge) ) uf = 1; 
00299   if ( !fit2.Fail() && (fit2.Rcs()<fit1.Rcs()) ) uf = 2;
00300 
00301   if (uf) {//two fit all ok
00302   //delete bad track
00303   MdcxHel fitme = (uf == 1) ? fit1 : fit2;
00304   }
00305   */
00306 }

void MdcMergeDups::dumpRecMdcTrack (   ) 
 

void MdcMergeDups::dumpRecMdcTrack (   ) 
  

00470                                   {
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     //std::cout<< "errmat   " << std::endl;
00507     //for (int i=0; i<15; i++){ std::cout<< " "<<tk->err(i); }
00508     //std::cout<< "   " << std::endl;
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     }//end of hit list
00519     std::cout << "  "<< std::endl;
00520   }//end of tk list
00521   std::cout << "  "<< std::endl;
00522 }

bool MdcMergeDups::eraseTdsTrack (  RecMdcTrackCol::iterator  tk  ) 
 

bool MdcMergeDups::eraseTdsTrack (  RecMdcTrackCol::iterator  tk  ) 
  

00455                                                          {
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     //hitList->erase(iterHit);
00464   }
00465   trackList->erase(tk);
00466   return true;
00467 }

StatusCode MdcMergeDups::execute (   ) 
 

StatusCode MdcMergeDups::execute (   ) 
  

00113                                  {
00114   MsgStream log(msgSvc(), name());
00115   log << MSG::INFO << "in execute()" << endreq;
00116   setFilterPassed(false); 
00117 
00118   m_bunchT0 = -999.;
00119   SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00120   if (!aevtimeCol || aevtimeCol->size()==0) {
00121     log << MSG::WARNING<< " Could not find RecEsTimeCol"<< endreq;
00122     return StatusCode::SUCCESS;
00123   }
00124 
00125   RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00126   for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00127     m_bunchT0 =  (*iter_evt)->getTest();
00128   }
00129 
00130 
00131   int nMerged = mergeCurl();
00132 
00133   if(m_debug>0) {
00134     std::cout<<name()<<": Merged "<<nMerged << " track "<<  std::endl;
00135     dumpRecMdcTrack();
00136   }
00137 
00138   return StatusCode::SUCCESS;   
00139 }

StatusCode MdcMergeDups::finalize (   ) 
 

StatusCode MdcMergeDups::finalize (   ) 
  

00142                                  {  
00143   MsgStream log(msgSvc(), name());
00144   log << MSG::INFO << "in finalize()" << endreq;        
00145 
00146   return StatusCode::SUCCESS;   
00147 }

StatusCode MdcMergeDups::initialize (   ) 
 

StatusCode MdcMergeDups::initialize (   ) 
  

00092                                    {  
00093   MsgStream log(msgSvc(), name());
00094   log << MSG::INFO << "in initialize()" << endreq;      
00095   StatusCode sc;
00096 
00097 
00098   //Initailize magnetic filed 
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   m_Bz = m_bfield->bFieldZ();
00107   log << MSG::INFO << "field z = "<<m_Bz<< endreq;
00108 
00109   return StatusCode::SUCCESS;   
00110 }

int MdcMergeDups::mergeCurl (  void   ) 
 

int MdcMergeDups::mergeCurl (  void   ) 
  

00151                            {
00152   //-------------------------------------------------------------------
00153 
00154   SmartDataPtr<RecMdcTrackCol> trackList(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00155   if (!trackList) return -1;
00156 
00157   int needMerge = 0;
00158 
00159   //...Merging. Search a track to be merged...
00160   RecMdcTrackCol::iterator iterRefTk = trackList->begin();
00161   for (; iterRefTk != trackList->end(); iterRefTk++) {
00162     RecMdcTrack* refTk = *iterRefTk;
00163     if (refTk->stat()<0) continue;
00164     std::vector<RecMdcTrack*> mergeTkList;
00165     mergeTkList.push_back(refTk);
00166 
00167 
00168     bool curl = false;
00169     int sameParm = 0;
00170     RecMdcTrackCol::iterator iterTestTk = trackList->begin();
00171     for (; iterTestTk != trackList->end(); iterTestTk++) {
00172       RecMdcTrack* testTk = *iterTestTk;
00173       if (iterRefTk == iterTestTk || (testTk->stat()<0)) continue;
00174 
00175       //-- overlapRatio cut 0.7 by jialk, original is 0.8
00176       if (testByOverlapHit(refTk,testTk)){
00177         if(m_debug>0)std::cout<<__FILE__<<" overlape tk:" <<refTk->trackId()<<" with "<<testTk->trackId()<<  std::endl;
00178         mergeTkList.push_back(testTk);
00179         curl = true;
00180       }
00181       sameParm = testByParam(refTk,testTk);
00182       if(sameParm >0) {
00183         if(m_debug>0) std::cout<<__FILE__<<" same param tk:" <<refTk->trackId()<<" with "<<testTk->trackId()<<  std::endl;
00184         mergeTkList.push_back(testTk);
00185       }
00186     } 
00187     if (mergeTkList.size()>1 && curl) needMerge = doMergeCurl(mergeTkList); 
00188     if ((needMerge < 999) && mergeTkList.size()>1 ) needMerge = doMergeLong(mergeTkList); 
00189     //if ((needMerge <999) && mergeTkList.size()==2 && (sameParm==2) ) needMerge = doMergeOdd(mergeTkList); //FIXME
00190   }
00191 
00192   //return 0 if No track need merged
00193   if( needMerge <=0 ) return 0;
00194 
00195   // reset track Id
00196   iterRefTk = trackList->begin();
00197   int iTk=0;
00198   int nDeleted = 0;
00199   for (; iterRefTk != trackList->end(); ) {
00200     if ( (*iterRefTk)->stat() >= 0 ){
00201       (*iterRefTk)->setTrackId(iTk);
00202       iterRefTk++;
00203       iTk++;
00204     }else {
00205       int id = (*iterRefTk)->trackId();
00206       bool erased = eraseTdsTrack(iterRefTk);
00207       if ( erased ){ 
00208         nDeleted++;
00209         if(m_debug>0)std::cout<<__FILE__<<" erase track No."<<id<<  std::endl;
00210       }else {
00211         if(m_debug>0)std::cout<<__FILE__<<" erase failed !"<<  std::endl;
00212       }
00213     }
00214 
00215   }
00216   if(m_debug>0) std::cout<<__FILE__<<" After merge save "<<iTk<<" tracks"<<  std::endl;
00217 
00218   return nDeleted;
00219 }

int MdcMergeDups::mergeDups (  void   ) 
 

int MdcMergeDups::mergeDups (  void   ) 
 

void MdcMergeDups::store (  TrkRecoTrk *  aTrack  ) 
 

void MdcMergeDups::store (  TrkRecoTrk *  aTrack  ) 
  

00439                                           {
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);//aTrack have been deleted in ~MdcTrack() 
00450   //tkStat: 0,Tsf 1,CurlFinder 2,PatRec 3,MdcxReco 4,MergeCurl
00451   int tkStat = 4;
00452   mdcTrack.storeTrack(-1, trackList, hitList, tkStat);
00453 }

int MdcMergeDups::testByOverlapHit (  RecMdcTrack *  refTk,
RecMdcTrack *  testTk
) 
 

int MdcMergeDups::testByOverlapHit (  RecMdcTrack *  refTk,
RecMdcTrack *  testTk
) 
  

00349                                                                          {
00350   //-------------------------------------------------------------------
00351   int overlaped = 0;
00352   if ((testTk->pxy() >= m_mergePt) || (refTk->pxy() >= m_mergePt)) return overlaped;
00353 
00354   HitRefVec testHits = testTk->getVecHits();
00355   int nHit = testHits.size();
00356   int nOverlap = 0;
00357 
00358   HitRefVec::iterator iterHit = testHits.begin();
00359   for (; iterHit != testHits.end(); iterHit++) {
00360     RecMdcHit* hit = *iterHit;
00361 
00362     //-- load for Axial and Stereo layer are 3,4 by jialk, original is 2,3
00363     double load = m_mergeLoadAx; 
00364     bool isStLayer = (m_gm->Layer(MdcID::layer(hit->getMdcId()))->view() != 0);
00365     if(isStLayer) load = m_mergeLoadSt;
00366 
00367     //helix parameters
00368     double vx0  = refTk->getVX0();
00369     double vy0  = refTk->getVY0();
00370     double vz0  = refTk->getVZ0();
00371     double dr   = refTk->helix(0);
00372     double phi0 = refTk->helix(1);
00373     double r    = 10000./ (Constants::c * m_Bz*refTk->helix(2));
00374     double dz   = refTk->helix(3);
00375     double tanl = refTk->helix(4);
00376 
00377     //center of circle
00378     double xc = vx0 + (dr + r) * cos(phi0);
00379     double yc = vy0 + (dr + r) * sin(phi0);
00380 
00381     //position of hit
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     //distance from center of circle to hit
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     //is this hit overlaped ? 
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 }

int MdcMergeDups::testByParam (  RecMdcTrack *  refTk,
RecMdcTrack *  testTk
) 
 

int MdcMergeDups::testByParam (  RecMdcTrack *  refTk,
RecMdcTrack *  testTk
) 
  

00224                                                                     {
00225   //-------------------------------------------------------------------
00226   int overlaped = 0;
00227 
00228 
00229   //Convert to Babar track convension
00230   double omega1 = (Constants::c * m_Bz*refTk->helix(2))/10000.;
00231   double omega2 = (Constants::c * m_Bz*testTk->helix(2))/10000.;
00232   //phi0_babar = phi0_belle + pi/2   [0,2pi)
00233   double phi01  = refTk->helix(1)+Constants::pi/2.;
00234   double phi02  = testTk->helix(1)+Constants::pi/2.;
00235   while(phi01>Constants::twoPi) phi01 -= Constants::twoPi;
00236   while(phi02>Constants::twoPi) phi02 -= Constants::twoPi;
00237   double d01    = -refTk->helix(0);
00238   double d02    = -testTk->helix(0);
00239   double dphi0  = fabs(phi01 - phi02);
00240   double dd0    = fabs(d01 - d02);
00241   double prodo  = omega1*omega2;
00242   double r1=100000.;
00243   double r2=100000.;
00244   if (fabs(omega1)>0.00001) r1 = 1.0/fabs(omega1);
00245   if (fabs(omega2)>0.00001) r2 = 1.0/fabs(omega2); //FIXME
00246   double pdrad = fabs((r1-r2)/(r1+r2)) ;
00247 
00248   if (2==m_debug){
00249     std::cout << "  fabs(d01 - d02)  " << fabs(d01 - d02) << std::endl;
00250     std::cout << "  fabs(phi01-phi02)  " << fabs(phi01-phi02) << std::endl;
00251   }
00252   //  Try to merge pair that looks like duplicates (same charge)
00253   if ( (prodo > 0.) && (dd0 < m_maxDd0InMerge) && (dphi0 < m_maxDphi0InMerge) &&
00254       (pdrad < m_maxPdradInMerge)) {
00255     overlaped = 1;
00256   }
00257 
00258   //  Try to merge pair that looks like albedo (opp charge, large d0)
00259   if ( (prodo < 0.) && (fabs(d01+d02) < 4.0) && (dd0 > 47.0) &&
00260       (fabs( dphi0 - Constants::pi) < m_maxDphi0InMerge) 
00261       && (pdrad < m_maxPdradInMerge)) {
00262     overlaped = 2;
00263   }
00264 
00265   return overlaped;
00266 }

Member Data Documentation

BField* MdcMergeDups::m_bfield [private]
 

BField* MdcMergeDups::m_bfield [private]
 

double MdcMergeDups::m_bunchT0 [private]
 

double MdcMergeDups::m_Bz [private]
 

int MdcMergeDups::m_debug [private]
 

const MdcDetector* MdcMergeDups::m_gm [private]
 

const MdcDetector* MdcMergeDups::m_gm [private]
 

double MdcMergeDups::m_maxDd0InMerge [private]
 

double MdcMergeDups::m_maxDphi0InMerge [private]
 

double MdcMergeDups::m_maxPdradInMerge [private]
 

double MdcMergeDups::m_maxRcsInMerge [private]
 

double MdcMergeDups::m_mergeLoadAx [private]
 

double MdcMergeDups::m_mergeLoadSt [private]
 

double MdcMergeDups::m_mergeOverlapRatio [private]
 

double MdcMergeDups::m_mergePt [private]
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 16:26:46 2011 for BOSS6.5.5 by  doxygen 1.3.9.1