MdcMergeDups Class Reference

#include <MdcMergeDups.h>

Inheritance diagram for MdcMergeDups:

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

Detailed Description

Definition at line 48 of file MdcMergeDups.h.

Constructor & Destructor Documentation

MdcMergeDups()

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

Definition at line 59 of file MdcMergeDups.cxx.

    : Algorithm( name, pSvcLocator ) {
  declareProperty( "debug", m_debug = 0 );
  // cuts for mergeDups()
  declareProperty( "maxDd0InMerge", m_maxDd0InMerge = 2.7 );
  declareProperty( "maxDphi0InMerge", m_maxDphi0InMerge = 0.15 );
  declareProperty( "maxDPdradInMerge", m_maxPdradInMerge = 0.22 );
  declareProperty( "maxRcsInMerge", m_maxRcsInMerge = 18. );
  // cuts for mergeCurl()
  declareProperty( "mergePt", m_mergePt = 0.13 );
  declareProperty( "mergeLoadAx", m_mergeLoadAx = 3. );
  declareProperty( "mergeLoadSt", m_mergeLoadSt = 4. );
  declareProperty( "mergeOverlapRatio", m_mergeOverlapRatio = 0.7 );
}

Referenced by MdcMergeDups().

~MdcMergeDups()

MdcMergeDups::~MdcMergeDups ( )

virtual

Definition at line 77 of file MdcMergeDups.cxx.

{ delete m_bfield; }

Member Function Documentation

beginRun()

StatusCode MdcMergeDups::beginRun

(

)

Definition at line 79 of file MdcMergeDups.cxx.

                                  {
  // Detector geometry
  m_gm = MdcDetector::instance();
  if ( NULL == m_gm ) return StatusCode::FAILURE;
  return StatusCode::SUCCESS;
}

Referenced by execute().

doMergeCurl()

int MdcMergeDups::doMergeCurl

(

std::vector< RecMdcTrack * >

mergeTkList

)

Definition at line 426 of file MdcMergeDups.cxx.

                                                                   {
  //-------------------------------------------------------------------
  int                                 innerMostTkId      = 999;
  RecMdcTrack*                        innerMostTk        = NULL;
  unsigned                            innerMostLayerOfTk = 999;
  std::vector<RecMdcTrack*>::iterator itTk               = mergeTkList.begin();
  for ( int iTk = 0; itTk != mergeTkList.end(); itTk++, iTk++ )
  {
    RecMdcTrack* tk             = ( *itTk );
    unsigned     innerMostLayer = 999;
    for ( unsigned iHit = 0; iHit < tk->getVecHits().size(); iHit++ )
    {
      unsigned layer = MdcID::layer( tk->getVecHits()[iHit]->getMdcId() );
      if ( layer < innerMostLayer ) innerMostLayer = layer;
    }
 
    if ( m_debug > 0 )
      std::cout << __FILE__ << " to be merged track id=" << tk->trackId() << std::endl;
    // test inner most layer id; if same, test dz
    if ( innerMostLayer < innerMostLayerOfTk )
    {
      innerMostTkId = iTk;
      innerMostTk   = tk;
    }
    else if ( innerMostLayer == innerMostLayerOfTk )
    {
      // test by dz
      if ( tk->helix( 3 ) < innerMostTk->helix( 3 ) )
      {
        innerMostTkId = iTk;
        innerMostTk   = tk;
      }
    }
  } // end of for mergeTkList
  innerMostTk->setStat( -1 );
 
  return innerMostTkId;
}

Referenced by mergeCurl().

doMergeLong()

int MdcMergeDups::doMergeLong

(

std::vector< RecMdcTrack * >

mergeTkList

)

Definition at line 287 of file MdcMergeDups.cxx.

                                                                   {
  //-------------------------------------------------------------------
  // merge hitlist
  double                              minRcs   = 999.;
  int                                 bestTkId = 999;
  RecMdcTrack*                        bestTk   = NULL;
  std::vector<RecMdcTrack*>::iterator itTk     = mergeTkList.begin();
  for ( int iTk = 0; itTk != mergeTkList.end(); itTk++, iTk++ )
  {
    RecMdcTrack* tk   = ( *itTk );
    double       chi2 = tk->chi2();
    double       ndf  = tk->ndof();
    if ( chi2 / ndf < minRcs )
    {
      bestTkId = tk->trackId();
      bestTk   = tk;
    }
  }
  if ( minRcs < m_maxRcsInMerge ) return bestTkId;
  bestTk->setStat( -1 );
 
  return 999;
  // FIXME
  /*
  //fit with track parameter respectively
  MdcxFittedHel fit1(dcxhlist, *iptr);
  MdcxFittedHel fit2(dcxhlist, *trkl[j]);
  int uf = 0;
  //get a best fit
  if ( !fit1.Fail() && (fit1.Rcs()<m_maxRcsInMerge) ) uf = 1;
  if ( !fit2.Fail() && (fit2.Rcs()<fit1.Rcs()) ) uf = 2;
 
  if (uf) {//two fit all ok
  //delete bad track
  MdcxHel fitme = (uf == 1) ? fit1 : fit2;
  }
  */
}

Referenced by mergeCurl().

dumpRecMdcTrack()

void MdcMergeDups::dumpRecMdcTrack

(

)

Definition at line 496 of file MdcMergeDups.cxx.

                                   {
  SmartDataPtr<RecMdcTrackCol> trackList( eventSvc(), EventModel::Recon::RecMdcTrackCol );
  if ( !trackList ) return;
  if ( trackList->size() != 4 ) setFilterPassed( true );
  std::cout << "N track after Merged = " << trackList->size() << std::endl;
  if ( m_debug <= 1 ) return;
  RecMdcTrackCol::iterator it = trackList->begin();
  for ( ; it != trackList->end(); it++ )
  {
    RecMdcTrack* tk = *it;
    std::cout << "//====RecMdcTrack " << tk->trackId() << "====:" << std::endl;
    cout << " d0 " << tk->helix( 0 ) << " phi0 " << tk->helix( 1 ) << " cpa " << tk->helix( 2 )
         << " z0 " << tk->helix( 3 ) << " tanl " << tk->helix( 4 ) << endl;
    std::cout << " q " << tk->charge() << " theta " << tk->theta() << " phi " << tk->phi()
              << " x0 " << tk->x() << " y0 " << tk->y() << " z0 " << tk->z() << " r0 "
              << tk->r() << endl;
    std::cout << " p " << tk->p() << " pt " << tk->pxy() << " px " << tk->px() << " py "
              << tk->py() << " pz " << tk->pz() << endl;
    std::cout << " tkStat " << tk->stat() << " chi2 " << tk->chi2() << " ndof " << tk->ndof()
              << " nhit " << tk->getNhits() << " nst " << tk->nster() << endl;
    // std::cout<< "errmat   " << std::endl;
    // for (int i=0; i<15; i++){ std::cout<< " "<<tk->err(i); }
    // std::cout<< "   " << std::endl;
 
    int nhits = tk->getVecHits().size();
    std::cout << nhits << " Hits: " << std::endl;
    for ( int ii = 0; ii < nhits; ii++ )
    {
      Identifier id( tk->getVecHits()[ii]->getMdcId() );
      int        layer = MdcID::layer( id );
      int        wire  = MdcID::wire( id );
      cout << "(" << layer << "," << wire << "," << tk->getVecHits()[ii]->getStat()
           << ",lr:" << tk->getVecHits()[ii]->getFlagLR() << ") ";
    } // end of hit list
    std::cout << "  " << std::endl;
  } // end of tk list
  std::cout << "  " << std::endl;
}

Referenced by execute().

eraseTdsTrack()

bool MdcMergeDups::eraseTdsTrack

(

RecMdcTrackCol::iterator

tk

)

Definition at line 481 of file MdcMergeDups.cxx.

                                                            {
  SmartDataPtr<RecMdcTrackCol> trackList( eventSvc(), EventModel::Recon::RecMdcTrackCol );
  if ( !trackList ) return false;
  SmartDataPtr<RecMdcHitCol> hitList( eventSvc(), EventModel::Recon::RecMdcHitCol );
  if ( !hitList ) return false;
  HitRefVec           hits    = ( *tk )->getVecHits();
  HitRefVec::iterator iterHit = hits.begin();
  for ( ; iterHit != hits.end(); iterHit++ )
  {
    // hitList->erase(iterHit);
  }
  trackList->erase( tk );
  return true;
}

Referenced by mergeCurl().

execute()

StatusCode MdcMergeDups::execute

(

)

Definition at line 107 of file MdcMergeDups.cxx.

                                 {
  if ( !m_beginRun )
  {
    StatusCode sc = beginRun();
    if ( sc.isFailure() )
    {
      error() << "beginRun failed" << endmsg;
      return StatusCode::FAILURE;
    }
    m_beginRun = true;
  }
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
  setFilterPassed( false );
 
  m_bunchT0 = -999.;
  SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
  if ( !aevtimeCol || aevtimeCol->size() == 0 )
  {
    log << MSG::WARNING << " Could not find RecEsTimeCol" << endmsg;
    return StatusCode::SUCCESS;
  }
 
  RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
  for ( ; iter_evt != aevtimeCol->end(); iter_evt++ ) { m_bunchT0 = ( *iter_evt )->getTest(); }
 
  int nMerged = mergeCurl();
 
  if ( m_debug > 0 )
  {
    std::cout << name() << ": Merged " << nMerged << " track " << std::endl;
    dumpRecMdcTrack();
  }
 
  return StatusCode::SUCCESS;
}

finalize()

StatusCode MdcMergeDups::finalize

(

)

Definition at line 145 of file MdcMergeDups.cxx.

                                  {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
 
  return StatusCode::SUCCESS;
}

initialize()

StatusCode MdcMergeDups::initialize

(

)

Definition at line 89 of file MdcMergeDups.cxx.

                                    {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
  StatusCode sc;
 
  // Initailize magnetic filed
  IBesMagFieldSvc* m_pIMF;
  sc = service( "MagneticFieldSvc", m_pIMF );
  if ( sc != StatusCode::SUCCESS )
  {
    log << MSG::ERROR << "Unable to open Magnetic field service" << endmsg;
    return StatusCode::FAILURE;
  }
  m_bfield = new BField( m_pIMF );
 
  return StatusCode::SUCCESS;
}

mergeCurl()

int MdcMergeDups::mergeCurl

(

void

)

Definition at line 153 of file MdcMergeDups.cxx.

                            {
  //-------------------------------------------------------------------
 
  SmartDataPtr<RecMdcTrackCol> trackList( eventSvc(), EventModel::Recon::RecMdcTrackCol );
  if ( !trackList ) return -1;
 
  int needMerge = 0;
 
  //...Merging. Search a track to be merged...
  RecMdcTrackCol::iterator iterRefTk = trackList->begin();
  for ( ; iterRefTk != trackList->end(); iterRefTk++ )
  {
    RecMdcTrack* refTk = *iterRefTk;
    if ( refTk->stat() < 0 ) continue;
    std::vector<RecMdcTrack*> mergeTkList;
    mergeTkList.push_back( refTk );
 
    bool                     curl       = false;
    int                      sameParm   = 0;
    RecMdcTrackCol::iterator iterTestTk = trackList->begin();
    for ( ; iterTestTk != trackList->end(); iterTestTk++ )
    {
      RecMdcTrack* testTk = *iterTestTk;
      if ( iterRefTk == iterTestTk || ( testTk->stat() < 0 ) ) continue;
 
      //-- overlapRatio cut 0.7 by jialk, original is 0.8
      if ( testByOverlapHit( refTk, testTk ) )
      {
        if ( m_debug > 0 )
          std::cout << __FILE__ << " overlape tk:" << refTk->trackId() << " with "
                    << testTk->trackId() << std::endl;
        mergeTkList.push_back( testTk );
        curl = true;
      }
      sameParm = testByParam( refTk, testTk );
      if ( sameParm > 0 )
      {
        if ( m_debug > 0 )
          std::cout << __FILE__ << " same param tk:" << refTk->trackId() << " with "
                    << testTk->trackId() << std::endl;
        mergeTkList.push_back( testTk );
      }
    }
    if ( mergeTkList.size() > 1 && curl ) needMerge = doMergeCurl( mergeTkList );
    if ( ( needMerge < 999 ) && mergeTkList.size() > 1 )
      needMerge = doMergeLong( mergeTkList );
    // if ((needMerge <999) && mergeTkList.size()==2 && (sameParm==2) ) needMerge =
    // doMergeOdd(mergeTkList); //FIXME
  }
 
  // return 0 if No track need merged
  if ( needMerge <= 0 ) return 0;
 
  // reset track Id
  iterRefTk    = trackList->begin();
  int iTk      = 0;
  int nDeleted = 0;
  for ( ; iterRefTk != trackList->end(); )
  {
    if ( ( *iterRefTk )->stat() >= 0 )
    {
      ( *iterRefTk )->setTrackId( iTk );
      iterRefTk++;
      iTk++;
    }
    else
    {
      int  id     = ( *iterRefTk )->trackId();
      bool erased = eraseTdsTrack( iterRefTk );
      if ( erased )
      {
        nDeleted++;
        if ( m_debug > 0 ) std::cout << __FILE__ << " erase track No." << id << std::endl;
      }
      else
      {
        if ( m_debug > 0 ) std::cout << __FILE__ << " erase failed !" << std::endl;
      }
    }
  }
  if ( m_debug > 0 )
    std::cout << __FILE__ << " After merge save " << iTk << " tracks" << std::endl;
 
  return nDeleted;
}

Referenced by execute().

mergeDups()

int MdcMergeDups::mergeDups

(

void

)

store()

void MdcMergeDups::store

(

TrkRecoTrk *

aTrack

)

Definition at line 465 of file MdcMergeDups.cxx.

                                             {
  SmartDataPtr<RecMdcTrackCol> trackList( eventSvc(), EventModel::Recon::RecMdcTrackCol );
  if ( !trackList ) return;
  SmartDataPtr<RecMdcHitCol> hitList( eventSvc(), EventModel::Recon::RecMdcHitCol );
  if ( !hitList ) return;
 
  assert( aTrack != NULL );
  TrkExchangePar helix = aTrack->fitResult()->helix( 0. );
 
  if ( m_debug > 1 ) std::cout << __FILE__ << " STORED" << std::endl;
  MdcTrack mdcTrack( aTrack ); // aTrack have been deleted in ~MdcTrack()
  // tkStat: 0,Tsf 1,CurlFinder 2,PatRec 3,MdcxReco 4,MergeCurl
  int tkStat = 4;
  mdcTrack.storeTrack( -1, trackList, hitList, tkStat );
}

testByOverlapHit()

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

Definition at line 366 of file MdcMergeDups.cxx.

                                                                            {
  //-------------------------------------------------------------------
  int overlaped = 0;
  if ( ( testTk->pxy() >= m_mergePt ) || ( refTk->pxy() >= m_mergePt ) ) return overlaped;
 
  HitRefVec testHits = testTk->getVecHits();
  int       nHit     = testHits.size();
  int       nOverlap = 0;
 
  HitRefVec::iterator iterHit = testHits.begin();
  for ( ; iterHit != testHits.end(); iterHit++ )
  {
    RecMdcHit* hit = *iterHit;
 
    //-- load for Axial and Stereo layer are 3,4 by jialk, original is 2,3
    double load      = m_mergeLoadAx;
    bool   isStLayer = ( m_gm->Layer( MdcID::layer( hit->getMdcId() ) )->view() != 0 );
    if ( isStLayer ) load = m_mergeLoadSt;
 
    // helix parameters
    double vx0  = refTk->getVX0();
    double vy0  = refTk->getVY0();
    double vz0  = refTk->getVZ0();
    double dr   = refTk->helix( 0 );
    double phi0 = refTk->helix( 1 );
    double Bz   = m_bfield->bFieldZ();
    double r    = 10000. / ( Constants::c * Bz * refTk->helix( 2 ) );
    double dz   = refTk->helix( 3 );
    double tanl = refTk->helix( 4 );
 
    // center of circle
    double xc = vx0 + ( dr + r ) * cos( phi0 );
    double yc = vy0 + ( dr + r ) * sin( phi0 );
 
    // position of hit
    double zHit = hit->getZhit();
    double phi  = ( vz0 + dz - zHit ) / ( r * tanl );
    double xHit = vx0 + dr * cos( phi0 ) + r * ( cos( phi0 ) - cos( phi0 + phi ) );
    double yHit = vy0 + dr * sin( phi0 ) + r * ( sin( phi0 ) - sin( phi0 + phi ) );
 
    // distance from center of circle to hit
    double dx          = xc - xHit;
    double dy          = yc - yHit;
    double dHit2Center = sqrt( dx * dx + dy * dy );
    double rTk         = fabs( r );
 
    // is this hit overlaped ?
    if ( ( dHit2Center > ( rTk - load ) ) && ( dHit2Center < ( rTk + load ) ) ) nOverlap++;
  }
 
  if ( nOverlap <= 0 ) return overlaped;
 
  double overlapRatio = double( nOverlap ) / double( nHit );
 
  if ( overlapRatio > m_mergeOverlapRatio ) overlaped = 1;
 
  return overlaped;
}

Referenced by mergeCurl().

testByParam()

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

Definition at line 242 of file MdcMergeDups.cxx.

                                                                       {
  //-------------------------------------------------------------------
  int overlaped = 0;
 
  // Convert to Babar track convension
  double Bz     = m_bfield->bFieldZ();
  double omega1 = ( Constants::c * Bz * refTk->helix( 2 ) ) / 10000.;
  double omega2 = ( Constants::c * Bz * testTk->helix( 2 ) ) / 10000.;
  // phi0_babar = phi0_belle + pi/2   [0,2pi)
  double phi01 = refTk->helix( 1 ) + Constants::pi / 2.;
  double phi02 = testTk->helix( 1 ) + Constants::pi / 2.;
  while ( phi01 > Constants::twoPi ) phi01 -= Constants::twoPi;
  while ( phi02 > Constants::twoPi ) phi02 -= Constants::twoPi;
  double d01   = -refTk->helix( 0 );
  double d02   = -testTk->helix( 0 );
  double dphi0 = fabs( phi01 - phi02 );
  double dd0   = fabs( d01 - d02 );
  double prodo = omega1 * omega2;
  double r1    = 100000.;
  double r2    = 100000.;
  if ( fabs( omega1 ) > 0.00001 ) r1 = 1.0 / fabs( omega1 );
  if ( fabs( omega2 ) > 0.00001 ) r2 = 1.0 / fabs( omega2 ); // FIXME
  double pdrad = fabs( ( r1 - r2 ) / ( r1 + r2 ) );
 
  if ( 2 == m_debug )
  {
    std::cout << "  fabs(d01 - d02)  " << fabs( d01 - d02 ) << std::endl;
    std::cout << "  fabs(phi01-phi02)  " << fabs( phi01 - phi02 ) << std::endl;
  }
  //  Try to merge pair that looks like duplicates (same charge)
  if ( ( prodo > 0. ) && ( dd0 < m_maxDd0InMerge ) && ( dphi0 < m_maxDphi0InMerge ) &&
       ( pdrad < m_maxPdradInMerge ) )
  { overlaped = 1; }
 
  //  Try to merge pair that looks like albedo (opp charge, large d0)
  if ( ( prodo < 0. ) && ( fabs( d01 + d02 ) < 4.0 ) && ( dd0 > 47.0 ) &&
       ( fabs( dphi0 - Constants::pi ) < m_maxDphi0InMerge ) && ( pdrad < m_maxPdradInMerge ) )
  { overlaped = 2; }
 
  return overlaped;
}

Referenced by mergeCurl().

---

The documentation for this class was generated from the following files:

• MdcMergeDups.h
• MdcMergeDups.cxx