FTFinder.cxx

Go to the documentation of this file.

#include <cmath>
#include <iostream>
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataManagerSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/IPartPropSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/SmartRef.h"
#include "GaudiKernel/SmartRefVector.h"
#include "GaudiKernel/StatusCode.h"
 
#include "EvTimeEvent/RecEsTime.h"
#include "EventModel/EventHeader.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "MdcGeomSvc/IMdcGeomSvc.h"
#include "MdcGeomSvc/MdcGeoLayer.h"
#include "MdcGeomSvc/MdcGeoSuper.h"
#include "MdcGeomSvc/MdcGeoWire.h"
#include "MdcRawEvent/MdcDigi.h"
#include "MdcRecEvent/RecMdcHit.h"
#include "MdcRecEvent/RecMdcTrack.h"
#include "RawDataProviderSvc/IRawDataProviderSvc.h"
#include "RawEvent/RawDataUtil.h"
#include "ReconEvent/ReconEvent.h"
#include "TrigEvent/TrigData.h"
#include "TrigEvent/TrigGTD.h"
#include "TrigEvent/TrigGTDProvider.h"
#include "TrigEvent/TrigMdc.h"
 
#include "FTFinder.h"
#include "FTGeom.h"
#include "FTLayer.h"
#include "FTList.h"
#include "FTSegment.h"
#include "FTSuperLayer.h"
#include "FTTrack.h"
#include "FTWire.h"
#include "MdcParameter.h"
#include "ntupleItem.h"
 
#ifndef OnlineMode
#  include "McTruth/McParticle.h"
 
#  include "AIDA/IAxis.h"
#  include "AIDA/IHistogram1D.h"
#  include "AIDA/IHistogram2D.h"
#  include "AIDA/IHistogram3D.h"
#  include "AIDA/IHistogramFactory.h"
#endif
 
using namespace Event;
 
#ifndef OnlineMode
 
extern NTuple::Item<long>   g_ntrkMC, g_eventNo;
extern NTuple::Array<float> g_theta0MC, g_phi0MC;
extern NTuple::Array<float> g_pxMC, g_pyMC, g_pzMC, g_ptMC;
 
// recon
 
extern NTuple::Item<long>   g_ntrk;
extern NTuple::Array<float> g_px, g_py, g_pz, g_pt, g_p;
extern NTuple::Array<float> g_phi, g_theta;
extern NTuple::Array<float> g_vx, g_vy, g_vz;
extern NTuple::Array<float> g_dr, g_phi0, g_kappa, g_dz, g_tanl;
extern NTuple::Item<float>  g_estime;
 
extern IHistogram1D* g_ncellMC;
extern IHistogram1D* g_ncell;
extern IHistogram1D* g_naxialhit;
extern IHistogram1D* g_nstereohit;
extern IHistogram1D* g_nhit;
extern IHistogram2D* g_hitmap[20];
 
extern int num_2Dtrk, num_3Dtrk, num_finaltrk;
 
#endif
 
FTSuperLayer*     FTFinder::superLayer( int id ) const { return m_geom->getSuperLayer( id ); }
FTList<FTTrack*>& FTFinder::tracks() { return _tracks; }
void              FTFinder::setAlgorithmPointer( Algorithm* alg ) { m_algorithm = alg; }
 
int FTFinder::getWireId( FTWire* w ) const { return w->getWireId(); }
 
float FTFinder::x2t( const FTLayer* l, const float x ) const {
  return ( x * x ) / ( param->_xtCoEff * param->_xtCoEff * l->csize() );
}
 
float FTFinder::t2x( const FTLayer* l, const float t ) const {
  float x = ( t > 0.0f ) ? (( param->_xtCoEff ))*sqrt( t * l->csize() ) : 0.0f;
  if ( x > 0.47f * l->csize() )
  {
    x = 0.0004f * t +
        0.47f * l->csize() *
            ( 1.0f - 0.0004f * 0.47f / ( ( param->_xtCoEff ) * ( param->_xtCoEff ) ) );
  }
 
  return x;
}
 
Hep3Vector FTFinder::vertex() const { return Hep3Vector( _vx, _vy, _vz ); }
 
MdcParameter* FTFinder::param = MdcParameter::instance();
 
FTFinder::FTFinder()
    : tOffSet( 0 )
    , t0Estime( -999 )
    , tEstFlag( 0 )
    , tEstime( 10 )
    , _tracks( 0 )
    , _linkedSegments( 0 )
    , _axialSegCollect( 0 )
    , _vx( 0 )
    , _vy( 0 )
    , _vz( 0 )
    , _ExpNo( 0 )
    , _RunNo( 0 )
    , m_total_trk( 0 )
    , pivot( 0, 0, 0 ) {
  StatusCode scmgn = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
  if ( scmgn != StatusCode::SUCCESS )
  { std::cout << "Unable to open Magnetic field service" << std::endl; }
}
 
void FTFinder::init() {
  if ( !param->_doIt ) return;
 
  // Get the Particle Properties Service
  IPartPropSvc* p_PartPropSvc;
  StatusCode    PartPropStatus = Gaudi::svcLocator()->service( "PartPropSvc", p_PartPropSvc );
 
  if ( !PartPropStatus.isSuccess() )
  {
    std::cerr << "Could not initialize Particle Properties Service" << std::endl;
    return;
  }
  m_particleTable = p_PartPropSvc->PDT();
 
  StatusCode RawData =
      Gaudi::svcLocator()->service( "RawDataProviderSvc", m_rawDataProviderSvc );
  if ( !RawData.isSuccess() )
  {
    std::cerr << "Could not load RawDataProviderSvc!" << m_rawDataProviderSvc << endmsg;
    return;
  }
}
 
void FTFinder::begin_run() {
  eventNo = 0;
  runNo   = 0;
  expflag = 0;
  if ( !param->_doIt ) return;
  _ExpNo = 0;
  _RunNo = 0;
 
  IMdcGeomSvc* mdcGeomSvc;
  StatusCode   sc = Gaudi::svcLocator()->service( "MdcGeomSvc", mdcGeomSvc );
 
  if ( !sc.isSuccess() ) return;
 
  m_geom = FTGeom::instance();
}
 
void FTFinder::event() {
  if ( !param->_doIt ) return;
 
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
 
  MsgStream log( msgSvc, "FTFinder" );
 
  IDataProviderSvc* eventSvc;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc ).ignore();
 
  //--
  // clear old information
  //--
  clear();
 
  // check whether  Recon  already registered
  DataObject* aReconEvent;
  eventSvc->findObject( "/Event/Recon", aReconEvent ).ignore();
  if ( !aReconEvent )
  {
    // register ReconEvent Data Object to TDS;
    ReconEvent* recevt = new ReconEvent;
    StatusCode  sc     = eventSvc->registerObject( "/Event/Recon", recevt );
    if ( sc != StatusCode::SUCCESS )
    {
      log << MSG::FATAL << "Could not register ReconEvent" << endmsg;
      return;
    }
  }
 
  // register Event start time
  IDataManagerSvc* dataManSvc;
  Gaudi::svcLocator()->service( "EventDataSvc", dataManSvc ).ignore();
 
  DataObject* aEsTimeEvent;
  eventSvc->findObject( "/Event/Recon/RecEsTimeCol", aEsTimeEvent ).ignore();
  if ( aEsTimeEvent )
  {
    dataManSvc->clearSubTree( "/Event/Recon/RecEsTimeCol" ).ignore();
    eventSvc->unregisterObject( "/Event/Recon/RecEsTimeCol" ).ignore();
  }
 
  RecEsTimeCol* aRecEsTimeCol = new RecEsTimeCol;
  StatusCode    est;
  est = eventSvc->registerObject( "/Event/Recon/RecEsTimeCol", aRecEsTimeCol );
  if ( est.isFailure() )
  {
    log << MSG::FATAL << "Could not register RecEsTimeCol" << endmsg;
    return;
  }
  log << MSG::DEBUG << "RecEsTimeCol registered successfully!" << endmsg;
 
  //--
  // update wirehit information
  //--
  updateMdc();
 
  //--
  // segment finding
  //--
  for ( int i = 0; i < 11; i++ ) { m_geom->getSuperLayer( i )->mkSegmentList(); }
 
  //--
  // reduce noise and get start time from segment fit
  //--
  for ( int i = 0; i < 10; i++ ) { m_geom->getSuperLayer( i )->reduce_noise( tEstime ); }
 
  getTestime();
 
  //--
  // axial segment linking
  //--
  mkTrackList();
 
  //--
  // 2D track fitting
  //--
  //(*_superLayer).reAppendSalvage();
  log << MSG::DEBUG << "number of 2D tracks: " << _tracks.size() << endmsg;
 
#ifndef OnlineMode
  num_2Dtrk += _tracks.size();
#endif
 
  for ( auto it = _tracks.begin(); it != _tracks.end(); )
  {
    if ( !( *it )->r_phiFit() )
    {
      delete *it;
      it            = _tracks.erase( it );
      int tmp_index = std::distance( _tracks.begin(), it );
      log << MSG::DEBUG << "===========> deleted 2D track :  " << tmp_index << endmsg;
    }
    else it++;
  }
 
  if ( t0Estime != -999 )
  {
    // begin to make Event start time
    RecEsTime* arecestime = new RecEsTime;
    arecestime->setTest( t0Estime );
    arecestime->setStat( tEstFlag );
    aRecEsTimeCol->push_back( arecestime );
  }
 
  if ( !_tracks.size() )
  {
    makeTds().ignore();
    return;
  }
 
  ///--
  /// vertex(r-phi) fit and event timing correction
  //--
  int vtx_flag = VertexFit( 0 );
  evtTiming    = ( param->_evtTimeCorr ) ? CorrectEvtTiming() : 0;
 
  ///--
  /// 2D track re-fitting
  //--
  if ( param->_hitscut == 1 )
  {
    for ( auto trk_i : _tracks ) { trk_i->r_phiReFit( _vx, _vy, vtx_flag ); }
  }
 
  //--
  // cut bad hits from 2D track re-fitting
  if ( param->_hitscut == 2 )
  {
    for ( auto trk_i : _tracks )
    {
      for ( int j = 0; j < 2; j++ )
      {
        i_rPhiFit = trk_i->r_phi2Fit( _vx, _vy, vtx_flag );
        if ( i_rPhiFit != -99 ) trk_i->r_phi3Fit( i_rPhiFit, _vx, _vy, vtx_flag );
      }
      trk_i->r_phi4Fit( _vx, _vy, vtx_flag );
    }
  }
 
  //--
  // stereo segment linking
  //--
  mkTrack3D();
 
  //--
  // final track fittng
  //--
  log << MSG::DEBUG << "number of 3D tracks: " << _tracks.size() << endmsg;
 
#ifndef OnlineMode
  num_3Dtrk += _tracks.size();
#endif
 
  for ( auto it = _tracks.begin(); it != _tracks.end(); )
  {
 
    int trk_index = std::distance( _tracks.begin(), it ) + 1;
 
#ifndef OnlineMode
    log << MSG::DEBUG << "=======> 3D track: " << trk_index << endmsg;
    ( *it )->printout();
#endif
 
    if ( ( *it )->get_nhits() < 18 )
    {
      delete ( *it );
      log << MSG::DEBUG << "================> deleted 3D track :  " << trk_index << endmsg;
      it = _tracks.erase( it );
    }
    else it++;
  }
 
  if ( !_tracks.size() )
  {
    makeTds().ignore();
    return;
  }
 
  for ( auto trk_i : _tracks ) trk_i->s_zFit();
 
  //--
  // find primary event vertex
  //--
  if ( param->_findEventVertex ) { VertexFit( 1 ); }
 
  log << MSG::DEBUG << "final number of tracks: " << _tracks.size() << endmsg;
 
#ifndef OnlineMode
  num_finaltrk += _tracks.size();
  if ( param->_mkMdst ) makeMdst();
#endif
 
  //--
  // output tracking result into TDS
  // by wangdy
  //--
  if ( param->_mkTds ) makeTds().ignore();
}
 
int FTFinder::updateMdc() {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
 
  MsgStream log( msgSvc, "FTFinder" );
 
  IDataProviderSvc* eventSvc;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc ).ignore();
 
  if ( !eventSvc )
  {
    log << MSG::FATAL << "Could not find EventDataSvc" << endmsg;
    return 0;
  }
 
  // Part 1: Get the event header, print out event and run number
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc, "/Event/EventHeader" );
  if ( !eventHeader )
  {
    log << MSG::FATAL << "Could not find Event Header" << endmsg;
    return 0;
  }
 
  log << MSG::DEBUG << "MdcFastTrkAlg: retrieved event: " << eventHeader->eventNumber()
      << "  run: " << eventHeader->runNumber() << endmsg;
 
  eventNo = eventHeader->eventNumber();
  runNo   = eventHeader->runNumber();
  if ( runNo > 0 ) expflag = 1;
  else expflag = 0;
  int digiId;
 
#ifndef OnlineMode
  g_eventNo = eventHeader->eventNumber();
 
  // Part 2: Retrieve MC truth
  SmartDataPtr<McParticleCol> mcParticleCol( eventSvc, "/Event/MC/McParticleCol" );
  if ( !mcParticleCol ) { log << MSG::WARNING << "Could not find McParticle" << endmsg; }
  else
  {
    McParticleCol::iterator iter_mc = mcParticleCol->begin();
 
    digiId     = 0;
    int ntrkMC = 0;
    int charge = 0;
 
    for ( ; iter_mc != mcParticleCol->end(); iter_mc++, digiId++ )
    {
      log << MSG::DEBUG << "MDC digit No: " << digiId << endmsg;
      log << MSG::DEBUG << " particleId = " << ( *iter_mc )->particleProperty() << endmsg;
      int statusFlags = ( *iter_mc )->statusFlags();
      int pid         = ( *iter_mc )->particleProperty();
      if ( pid > 0 )
      {
        if ( m_particleTable->particle( pid ) )
        { charge = (int)m_particleTable->particle( pid )->charge(); }
      }
      else if ( pid < 0 )
      {
        if ( m_particleTable->particle( -pid ) )
        {
          charge = (int)m_particleTable->particle( -pid )->charge();
          charge *= -1;
        }
      }
      else { log << MSG::WARNING << "wrong particle id, please check data" << endmsg; }
 
      if ( charge == 0 || abs( cos( ( *iter_mc )->initialFourMomentum().theta() ) ) > 0.93 )
        continue;
 
      g_theta0MC[ntrkMC] = ( *iter_mc )->initialFourMomentum().theta();
      g_phi0MC[ntrkMC]   = ( *iter_mc )->initialFourMomentum().phi();
      g_pxMC[ntrkMC]     = ( *iter_mc )->initialFourMomentum().px();
      g_pyMC[ntrkMC]     = ( *iter_mc )->initialFourMomentum().py();
      g_pzMC[ntrkMC]     = ( *iter_mc )->initialFourMomentum().pz();
      g_ptMC[ntrkMC] =
          0.001 * sqrt( g_pxMC[ntrkMC] * g_pxMC[ntrkMC] + g_pyMC[ntrkMC] * g_pyMC[ntrkMC] );
      ntrkMC++;
    }
    g_ntrkMC = ntrkMC;
  }
#endif
 
#ifdef MultiThread
  // Part 3: Retrieve MDC digi
  SmartDataPtr<MdcDigiCol> mdcDigiVec( eventSvc, "/Event/Digi/MdcDigiCol" );
  if ( !mdcDigiVec )
  {
    log << MSG::FATAL << "Could not find MDC digi" << endmsg;
    return ( StatusCode::FAILURE );
  }
  MdcDigiCol::iterator iter1 = mdcDigiVec->begin();
  digiId                     = 0;
  Identifier mdcId;
  int        layerId;
  int        wireId;
  for ( ; iter1 != mdcDigiVec->end(); iter1++, digiId++ )
  {
#else
  //--use RawDataProvider to get MdcDigi
  //  bool m_dropRecHits=true;//or false
  MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(); // call MdcRawDataProvider
  MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
  digiId                     = 0;
  Identifier mdcId;
  int        layerId;
  int        wireId;
  for ( ; iter1 != mdcDigiVec.end(); iter1++, digiId++ )
  {
#endif
 
    mdcId   = ( *iter1 )->identify();
    layerId = MdcID::layer( mdcId );
    wireId  = MdcID::wire( mdcId );
 
#ifndef OnlineMode
    log << MSG::DEBUG << "MDC digit No: " << digiId << endmsg;
    log << MSG::DEBUG
        << " time_channel = " << RawDataUtil::MdcTime( ( *iter1 )->getTimeChannel() )
        << " charge_channel = " << RawDataUtil::MdcCharge( ( *iter1 )->getChargeChannel() )
        << " layerId = " << layerId << " wireId = " << wireId << endmsg;
#endif
 
    const int localwid = wireId;
    const int wid      = m_geom->layer2wireStart( layerId ) + localwid;
 
#ifndef OnlineMode
    g_ncellMC->fill( wid, 1.0 );
#endif
 
    //... skip invalid wireID's
    if ( wid < 0 || wid > 6795 )
    {
      std::cout << "FTFinder::updateMdc(): wid out of range: " << wid << std::endl;
      continue;
    }
 
    FTWire* w   = m_geom->getWire( wid );
    float   tdc = RawDataUtil::MdcTime( ( *iter1 )->getTimeChannel() );
    float   adc = RawDataUtil::MdcCharge( ( *iter1 )->getChargeChannel() );
 
#ifndef OnlineMode
    if ( eventNo == 466 ) { g_hitmap[0]->fill( w->x(), w->y() ); }
    if ( eventNo == 538 ) { g_hitmap[1]->fill( w->x(), w->y() ); }
    if ( eventNo == 408 ) { g_hitmap[2]->fill( w->x(), w->y() ); }
    if ( eventNo == 499 ) { g_hitmap[3]->fill( w->x(), w->y() ); }
    if ( eventNo == 603 ) { g_hitmap[4]->fill( w->x(), w->y() ); }
    if ( eventNo == 938 ) { g_hitmap[5]->fill( w->x(), w->y() ); }
#endif
 
    // tdc/adc calibration
    // new tdc include drift time and flight time
    float time = tdc - sqrt( w->x() * w->x() + w->y() * w->y() ) / 30;
    w->time( time );
    w->wireId( wid );
    w->setAdc( adc );
 
    //... check adc validation
    FTLayer* lyr = w->layer();
 
    //... save to the array
    // w->distance(time*40*0.0001);  //original
    w->distance( 0.0 ); // by X.-R. Lu
    w->setChi2( 999 );
    w->state( FTWireHit );
    FTSuperLayer* sup = m_geom->getSuperLayer( m_geom->layer2superLayer( lyr->layerId() ) );
    sup->appendHit( w );
  }
 
  return 1;
}
 
void FTFinder::clear() {
  evtTiming = 0;
 
  for ( int i = 0; i < 11; i++ ) { m_geom->getSuperLayer( i )->clear(); }
 
  // for ( auto trk_i : _tracks ) delete trk_i;
  _tracks.clear();
 
  // for ( auto seg_i : _axialSegCollect ) delete seg_i;
  _axialSegCollect.clear();
 
  _vx = -99999.;
  _vy = -99999.;
  _vz = -99999.;
 
  for ( auto& v : tEstime ) { v.clear(); }
}
 
void FTFinder::getTestime() {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
 
  MsgStream log( msgSvc, "FTFinder" );
 
  IDataProviderSvc* eventSvc;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc ).ignore();
 
  float sumT = 0, estime = 0;
  int   n  = 0;
  t0Estime = -999;
 
  for ( auto& tEstimeVec : tEstime )
  {
    for ( auto tEstime_i : tEstimeVec )
    {
      if ( tEstime_i != 0 )
      {
        sumT += tEstime_i;
        n++;
      }
    }
  }
 
  if ( n != 0 )
  {
    estime = sumT / n;
    estime += _t0cal;
 
    if ( estime > -1 )
    {
      if ( expflag == 0 )
      {
        int nbunch = ( (int)( estime - tOffSet ) ) / _bunchtime;
        if ( ( (int)( estime - tOffSet ) ) % (int)_bunchtime > _bunchtime / 2 )
          nbunch = nbunch + 1;
        t0Estime = nbunch * _bunchtime + tOffSet;
      }
      else { t0Estime = estime; }
 
      int trigtimming = 0;
 
      // Retrieve trigger timming information
      // timing system: TOF:1,  MDC:2, EMC:3, NONE:0
      SmartDataPtr<TrigData> trigData( eventSvc, "/Event/Trig/TrigData" );
      if ( trigData )
      {
        trigtimming = trigData->getTimingType();
        log << MSG::INFO << "Timing type: " << trigData->getTimingType() << endmsg;
      }
 
      if ( trigtimming == 1 ) tEstFlag = 117;
      if ( trigtimming == 2 ) tEstFlag = 127;
      if ( trigtimming == 3 ) tEstFlag = 137;
      if ( trigtimming == 0 ) tEstFlag = 107;
 
#ifndef OnlineMode
      g_estime = estime;
#endif
    }
  }
}
 
void FTFinder::mkTrackList() {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
 
  MsgStream log( msgSvc, "FTFinder" );
 
  FTList<FTSegment*> inner_segments;
  FTList<FTSegment*> outer_segments;
 
  linkAxialSuperLayer234( inner_segments );
  linkAxialSuperLayer910( outer_segments );
 
  _axialSegCollect.insert( _axialSegCollect.end(), inner_segments.begin(),
                           inner_segments.end() );
  _axialSegCollect.insert( _axialSegCollect.end(), outer_segments.begin(),
                           outer_segments.end() );
 
  int innerN = inner_segments.size();
  int outerN = outer_segments.size();
 
  log << MSG::DEBUG << innerN << " segments in inner axial layers!" << endmsg;
  log << MSG::DEBUG << outerN << " segments in outer axial layers!" << endmsg;
 
  for ( auto inner : inner_segments )
  {
    if ( inner->track() ) continue;
    FTTrack* track_cand = nullptr;
 
#ifndef OnlineMode
    inner->printout();
#endif
 
    if ( outerN == 0 )
    {
      track_cand = linkAxialSegments( inner, nullptr );
      if ( !track_cand ) continue;
      _linkedSegments = FTList<FTSegment*>( 5 );
 
      // append this to _tracks
      const FTList<FTSegment*>& segments = track_cand->axial_segments();
 
      for ( auto seg_i : segments ) { seg_i->track( track_cand ); }
 
      _tracks.push_back( track_cand );
      track_cand->setFTFinder( this );
      continue;
    }
 
    for ( auto outer : outer_segments )
    {
      if ( outer->track() )
      { // already used
        if ( outer == outer_segments.back() )
        {                                 // if the last outer segment
          if ( inner->track() ) continue; // inner segment never used
          track_cand = linkAxialSegments( inner, nullptr );
          if ( !track_cand ) continue;
          _linkedSegments = FTList<FTSegment*>( 5 );
          // append this to _tracks
          const FTList<FTSegment*>& segments = track_cand->axial_segments();
 
          for ( auto seg_i : segments ) { seg_i->track( track_cand ); }
          _tracks.push_back( track_cand );
          track_cand->setFTFinder( this );
        }
        continue;
      }
 
      track_cand = linkAxialSegments( inner, outer );
 
      if ( !track_cand )
      { // link failed
        if ( outer == outer_segments.back() )
        {                                 // if the last outer segment
          if ( inner->track() ) continue; // inner segment never used
          track_cand = linkAxialSegments( inner, nullptr );
        }
        if ( !track_cand ) continue;
      }
 
      _linkedSegments = FTList<FTSegment*>( 5 );
 
      // compair this to appended track candidate
 
      const FTList<FTSegment*>& segments = track_cand->axial_segments();
 
      // compare this to the appended track candidate
      if ( inner->track() )
      { // already used
        FTTrack* trk_tmp = inner->track();
 
        int  nTrks    = _tracks.size();
        auto cache_it = _tracks.end();
 
        for ( auto it = _tracks.begin(); it != _tracks.end(); it++ )
        {
          if ( *it == trk_tmp )
          {
            cache_it = it;
            break;
          }
        }
 
        const FTList<FTSegment*>& segments2 = trk_tmp->axial_segments();
 
        int n2 = segments2.size();
        if ( segments.size() >= segments2.size() &&
             track_cand->chi2_kappa_tmp() < trk_tmp->chi2_kappa_tmp() )
        {
          // swap this and appended one
          for ( auto seg2_i : segments2 ) { seg2_i->track( nullptr ); }
          for ( auto seg_i : segments ) { seg_i->track( track_cand ); }
 
          *cache_it = track_cand;
          track_cand->setFTFinder( this );
          delete trk_tmp;
        }
        else { delete track_cand; }
      }
      else
      {
        // append this
        for ( auto seg_i : segments ) { seg_i->track( track_cand ); }
        _tracks.push_back( track_cand );
        track_cand->setFTFinder( this );
      }
    } // end of loop of outerN
  } // end of loop innerN
 
#ifndef OnlineMode
  int i = 0;
  for ( auto trk_i : _tracks )
  {
    const FTList<FTSegment*>& segments = trk_i->axial_segments();
    log << MSG::DEBUG << " loop connected axial track: " << i << endmsg;
    int l = segments.size();
    for ( int j = 0; j < l; j++ ) { segments[j]->printout(); }
    i++;
  }
#endif
}
 
FTTrack* FTFinder::linkAxialSegments( FTSegment* inner, FTSegment* outer ) {
  float chi2_kappa = param->_chi2_kappa;
  _linkedSegments.clear();
 
  if ( !outer )
  {
    // allow only 2, 3, 4,  axial segments in one track
    int n = inner->wireHits().size();
    _linkedSegments.push_back( inner );
 
    if ( n >= 7 ) return ( new FTTrack( _linkedSegments, inner->kappa(), chi2_kappa ) );
    else return nullptr;
  }
 
  _linkedSegments.push_back( outer );
 
  int   n       = _linkedSegments.size();
  float SigmaK  = outer->kappa();
  float SigmaRR = outer->r();
  SigmaRR *= SigmaRR;
  float SigmaKRR  = SigmaK * SigmaRR;
  float SigmaKKRR = SigmaK * SigmaKRR;
 
  FTSegment* s            = _linkedSegments[n - 1];
  FTSegment* innerSegment = nullptr;
 
  float SigmaK_cache, SigmaRR_cache, SigmaKRR_cache, SigmaKKRR_cache;
 
  float Min_chi2    = param->_Min_chi2;
  float inX         = s->incomingX();
  float inY         = s->incomingY();
  float in_r        = s->innerBoundHits().front()->layer()->r();
  float incomingPhi = s->incomingPhi();
 
  FTSegment*    next              = inner;
  const FTWire* NextOuterBoundHit = next->outerBoundHits().front();
 
  float outX = next->outgoingX();
  float outY = next->outgoingY();
 
  //////////////////////////////////////////////////////
  float _trk_d    = -2 * ( -1. / 2.99792458 / m_pmgnIMF->getReferField() ) / s->kappa();
  float _angle1   = asin( NextOuterBoundHit->layer()->r() / _trk_d );
  float _angle2   = asin( s->outerBoundHits().front()->layer()->r() / _trk_d );
  float _ang_diff = _angle2 - _angle1;
  float _diff     = s->outgoingPhi() - next->outgoingPhi();
  _diff           = _diff - ( int( _diff / M_PI ) ) * 2 * M_PI;
  float _require  = _ang_diff - _diff;
 
  // cut of connecting inner and outer segments
  if ( _require < -0.10 || _require > 0.11 ) return nullptr;
 
  float SegK  = next->kappa();
  float SegRR = next->r();
  SegRR *= SegRR;
  const float out_r = NextOuterBoundHit->layer()->r();
  // kappa = -2. * alpha * ((Vout X Vin)_z / |Vin|*|Vout|) / |Vin-Vout|
  float GapK = 2. * ( -1. / 2.99792458 / m_pmgnIMF->getReferField() ) *
               ( inX * outY - inY * outX ) /
               ( in_r * out_r *
                 sqrt( ( inX - outX ) * ( inX - outX ) + ( inY - outY ) * ( inY - outY ) ) );
  // float GapRR = (currentLayer==j+1||currentLayer==j+2) ? 0.5*(in_r+out_r) : in_r+out_r;
  float GapRR = 0.5 * ( in_r + out_r );
  GapRR *= GapRR;
  float SigmaK_tmp    = ( SigmaK + SegK + GapK );
  float SigmaRR_tmp   = SigmaRR + SegRR + GapRR;
  float SigmaKRR_tmp  = SigmaKRR + SegK * SegRR + GapK * GapRR;
  float SigmaKKRR_tmp = SigmaKKRR + SegK * SegK * SegRR + GapK * GapK * GapRR;
  float MuK_tmp       = SigmaK_tmp / ( 2 * n + 1 );
  float chi2 =
      ( MuK_tmp * MuK_tmp * SigmaRR_tmp - 2. * MuK_tmp * SigmaKRR_tmp + SigmaKKRR_tmp ) /
      ( 2 * n + 1 );
  if ( ( chi2 - chi2_kappa ) < Min_chi2 )
  {
    Min_chi2        = chi2;
    innerSegment    = next;
    SigmaK_cache    = SigmaK_tmp;
    SigmaRR_cache   = SigmaRR_tmp;
    SigmaKRR_cache  = SigmaKRR_tmp;
    SigmaKKRR_cache = SigmaKKRR_tmp;
  }
  if ( innerSegment )
  {
    _linkedSegments.push_back( inner );
    n          = _linkedSegments.size();
    SigmaK     = SigmaK_cache;
    SigmaRR    = SigmaRR_cache;
    SigmaKRR   = SigmaKRR_cache;
    SigmaKKRR  = SigmaKKRR_cache;
    chi2_kappa = Min_chi2;
 
    if ( n > 1 ) return ( new FTTrack( _linkedSegments, SigmaK / ( 2 * n - 1 ), chi2_kappa ) );
  }
  else return nullptr;
  return nullptr;
}
 
void FTFinder::linkAxialSuperLayer234( FTList<FTSegment*>& inner_segments ) {
  FTList<FTSegment*>  _segments34;
  FTList<FTSegment*>& SuperLayer3Segments = m_geom->getSuperLayer( 3 )->segments();
  FTList<FTSegment*>& SuperLayer4Segments = m_geom->getSuperLayer( 4 )->segments();
 
  linkAxialSegments_step( SuperLayer3Segments, SuperLayer4Segments, _segments34,
                          param->_D_phi2, param->_chi2_2 );
  _segments34.insert( _segments34.end(), SuperLayer3Segments.begin(),
                      SuperLayer3Segments.end() );
  SuperLayer3Segments.clear();
 
  FTList<FTSegment*>& SuperLayer2Segments = m_geom->getSuperLayer( 2 )->segments();
  linkAxialSegments_step( SuperLayer2Segments, _segments34, inner_segments, param->_D_phi1,
                          param->_chi2_1 );
  inner_segments.insert( inner_segments.end(), _segments34.begin(), _segments34.end() );
 
  // zoujh: connect the 2&4 superlayers leap over 3
 
  for ( auto inner_it = SuperLayer2Segments.begin(); inner_it != SuperLayer2Segments.end(); )
  {
    FTSegment* inner = *inner_it;
 
    bool erased_inner = false; // a flag to indicate whether the inner segment is erased
 
    float in_outerPhi = inner->outgoingPhi();
    for ( auto outer_it = SuperLayer4Segments.begin(); outer_it != SuperLayer4Segments.end(); )
    {
      FTSegment* outer = *outer_it;
 
      float out_innerPhi = outer->incomingPhi();
      float D_phi        = fabs( in_outerPhi - out_innerPhi );
      if ( D_phi > M_PI ) D_phi = 2 * M_PI - D_phi;
 
      /////////////////////////////////////////////////////
 
      if ( D_phi > M_PI / 12.5 )
      {
        outer_it++;
        continue;
      }
 
      float inX   = inner->outgoingX();
      float inY   = inner->outgoingY();
      float outX  = outer->incomingX();
      float outY  = outer->incomingY();
      float in_r  = inner->outerBoundHits().front()->layer()->r();
      float out_r = outer->innerBoundHits().front()->layer()->r();
      float GapK =
          2. * ( -1. / 2.99792458 / m_pmgnIMF->getReferField() ) *
          ( inY * outX - inX * outY ) /
          ( in_r * out_r *
            sqrt( ( inX - outX ) * ( inX - outX ) + ( inY - outY ) * ( inY - outY ) ) );
      // float cache_in = ((inner->kappa()+outer->kappa()+GapK)/3.0 - inner->kappa())/in_r;
      float cache_in = ( ( outer->kappa() + GapK ) / 3.0 - inner->kappa() * 2.0 / 3.0 ) / in_r;
      float cache_out =
          ( ( inner->kappa() + GapK ) / 3.0 - outer->kappa() * 2.0 / 3.0 ) / out_r;
      float cache_gap = ( ( inner->kappa() + outer->kappa() ) / 3.0 - GapK * 2.0 / 3.0 ) *
                        2.0 / ( in_r + out_r );
      float chi2_z = cache_in * cache_in + cache_out * cache_out + cache_gap * cache_gap;
 
      if ( chi2_z > param->_D_phi1 )
      {
        outer_it++;
        continue;
      }
 
      /////////////////////////////////////////////////////
      inner->connect_outer( outer->wireHits(), outer->outerBoundHits() );
      inner->update();
      inner_segments.push_back( inner );
 
      inner_it     = SuperLayer2Segments.erase( inner_it );
      erased_inner = true;
 
      delete *outer_it;
      outer_it = SuperLayer4Segments.erase( outer_it );
 
      break;
    }
    if ( !erased_inner ) inner_it++;
  }
 
  inner_segments.insert( inner_segments.end(), SuperLayer2Segments.begin(),
                         SuperLayer2Segments.end() );
  inner_segments.insert( inner_segments.end(), SuperLayer4Segments.begin(),
                         SuperLayer4Segments.end() );
}
 
void FTFinder::linkAxialSuperLayer910( FTList<FTSegment*>& outer_segments ) {
  FTList<FTSegment*>& SuperLayer9Segments  = m_geom->getSuperLayer( 9 )->segments();
  FTList<FTSegment*>& SuperLayer10Segments = m_geom->getSuperLayer( 10 )->segments();
  linkAxialSegments_step( SuperLayer9Segments, SuperLayer10Segments, outer_segments,
                          param->_D_phi3, param->_chi2_3 );
 
  outer_segments.insert( outer_segments.end(), SuperLayer9Segments.begin(),
                         SuperLayer9Segments.end() );
 
  outer_segments.insert( outer_segments.end(), SuperLayer10Segments.begin(),
                         SuperLayer10Segments.end() );
}
 
void FTFinder::linkAxialSegments_step( FTList<FTSegment*>& InnerSegments,
                                       FTList<FTSegment*>& OuterSegments,
                                       FTList<FTSegment*>& Connected, float maxDphi,
                                       float maxChi2 ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
  MsgStream log( msgSvc, "FTFinder" );
 
  for ( auto inner_it = InnerSegments.begin();
        inner_it != InnerSegments.end(); ) // will erase element in loop, so no "inner_it++"
                                           // here
  {
    FTSegment* inner = *inner_it;
 
    const FTLayer* in_outerBound     = inner->outerBoundHits().front()->layer();
    float          in_outerPhi       = inner->outgoingPhi();
    float          min_Dphi          = M_PI / 2;
    auto           min_Dphi_outer_it = OuterSegments.end();
 
    for ( auto outer_it = OuterSegments.begin(); outer_it != OuterSegments.end(); outer_it++ )
    {
      FTSegment* outer = *outer_it;
      float      D_phi = fabs( in_outerPhi - outer->incomingPhi() );
      if ( D_phi > M_PI ) D_phi = 2 * M_PI - D_phi;
 
      if ( D_phi > min_Dphi ) continue;
 
      float inX   = inner->incomingX();
      float inY   = inner->incomingY();
      float outX  = outer->outgoingX();
      float outY  = outer->outgoingY();
      float in_r  = inner->innerBoundHits().front()->layer()->r();
      float out_r = outer->outerBoundHits().front()->layer()->r();
      float allK =
          2. * ( -1. / 2.99792458 / m_pmgnIMF->getReferField() ) *
          ( inY * outX - inX * outY ) /
          ( in_r * out_r *
            sqrt( ( inX - outX ) * ( inX - outX ) + ( inY - outY ) * ( inY - outY ) ) );
      float cache_in  = ( ( outer->kappa() + allK ) / 3.0 - inner->kappa() * 2 / 3.0 ) / in_r;
      float cache_out = ( ( inner->kappa() + allK ) / 3.0 - outer->kappa() * 2 / 3.0 ) / out_r;
      float cache_all = ( ( inner->kappa() + outer->kappa() ) / 3.0 - allK * 2 / 3.0 ) * 2.0 /
                        ( in_r + out_r );
      float chi2_z = cache_in * cache_in + cache_out * cache_out + cache_all * cache_all;
 
      log << MSG::DEBUG << "D_phi: " << D_phi << " chi2_z: " << chi2_z
          << " maxChi2: " << maxChi2 << endmsg;
 
      if ( chi2_z > maxChi2 ) continue;
      min_Dphi          = D_phi;
      min_Dphi_outer_it = outer_it;
    }
 
    if ( min_Dphi_outer_it == OuterSegments.end() )
    {
      inner_it++;
      continue;
    };
 
    log << MSG::DEBUG << "min_Dphi: " << min_Dphi << endmsg;
 
    FTSegment* outer = *min_Dphi_outer_it;
 
    const FTLayer* out_innerBound = outer->innerBoundHits().front()->layer();
 
    int dLayerId = out_innerBound->layerId() - in_outerBound->layerId();
 
    if ( dLayerId == 1 && min_Dphi > maxDphi )
    {
      inner_it++;
      continue;
    };
    if ( dLayerId == 2 && min_Dphi > maxDphi * 1.5 )
    {
      inner_it++;
      continue;
    };
    if ( dLayerId == 3 && min_Dphi > maxDphi * 2.25 )
    {
      inner_it++;
      continue;
    };
    if ( dLayerId != 1 && dLayerId != 2 && dLayerId != 3 )
    {
      inner_it++;
      continue;
    };
 
    inner->connect_outer( outer->wireHits(), outer->outerBoundHits() );
    inner->update();
    Connected.push_back( inner );
    inner_it = InnerSegments.erase( inner_it );
 
    delete outer;
    OuterSegments.erase( min_Dphi_outer_it );
 
    log << MSG::DEBUG << "DONE!!" << endmsg;
  }
}
 
void FTFinder::mkTrack3D() {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
  MsgStream log( msgSvc, "FTFinder" );
 
  FTList<FTSegment*> multi_trk_cand;
 
  for ( int i = 8; i >= 0; i-- )
  {
    if ( i == 4 ) i -= 3;
    FTList<FTSegment*>& segments = m_geom->getSuperLayer( i )->segments();
 
    for ( auto s : segments )
    {
 
#ifndef OnlineMode
      s->printout();
#endif
 
      int      nTrack = 0;
      FTTrack* t;
 
      int k = 0;
      for ( auto tmp_trk : _tracks )
      {
#ifndef OnlineMode
        log << MSG::DEBUG << "coupling to track: " << k << endmsg;
#endif
        if ( s->update3D( tmp_trk ) )
        {
          // calculate s and z
          t = tmp_trk;
          nTrack++;
        }
        k++;
      }
 
      switch ( nTrack )
      {
      case 0: continue;
      case 1: t->append_stereo_cache( s ); break;
      default: multi_trk_cand.push_back( s ); break;
      }
    }
    // whose relation between this is unique
    for ( auto t : _tracks ) t->updateSZ();
  }
 
  // link segments by tanLambda & dz
  for ( auto t : _tracks ) t->linkStereoSegments();
  for ( auto t : multi_trk_cand ) t->linkStereoSegments();
}
 
int FTFinder::VertexFit2D() {
  int nTrks = _tracks.size();
  if ( nTrks < 2 )
  {
    _vx = -99999.;
    _vy = -99999.;
    _vz = -99999.;
    return 0;
  }
 
  FTList<float>  px;
  FTList<float>  py;
  FTList<float>  dx;
  FTList<float>  dy;
  FTList<double> weight;
 
  for ( auto t : _tracks )
  {
    const Lpav& la = t->lpav();
    HepVector   a  = la.Hpar( pivot );
 
    const float dr_i = a( 1 );
    if ( fabs( a( 1 ) ) > 1.5 ) continue;
    const float px_i = -sin( a( 2 ) );
    const float py_i = cos( a( 2 ) );
 
    float weight_i = la.chisq() / ( la.nc() * 0.02 );
 
    px.push_back( px_i );
    py.push_back( py_i );
    dx.push_back( dr_i * py_i );
    dy.push_back( -dr_i * px_i );
    weight.push_back( exp( -weight_i * weight_i ) );
  }
 
  if ( dx.size() < 2 )
  {
    _vx = -99999.;
    _vy = -99999.;
    _vz = -99999.;
    return 0;
  }
 
  double S_weight = 0.;
  for ( int i = dx.size() - 1; i; i-- )
  {
    const float px_i     = px[i];
    const float py_i     = py[i];
    const float dx_i     = dx[i];
    const float dy_i     = dy[i];
    const float weight_i = weight[i];
    for ( int j = 0; j < i; j++ )
    {
      const float px_j = px[j];
      const float py_j = py[j];
 
      const float ddx = dx[j] - dx_i;
      const float ddy = dx[j] - dy_i;
 
      const float tmp_par = px_i * py_j - px_j * py_i;
      const float par     = ( py_j * ddx - px_j * ddy ) / tmp_par;
 
      double weight_ij = weight_i * weight[j];
      S_weight += weight_i * weight[j];
      if ( tmp_par == 0 || par < -0.5 || ( py_i * ddx - px_i * ddy ) / tmp_par < -0.5 ||
           fabs( ( px_i * px_j + py_i * py_j ) /
                 sqrt( ( px_i * px_i + py_i * py_i ) * ( px_j * px_j + py_j * py_j ) ) ) >
               0.86 )
      {
        _vx += ( dx_i + 0.5 * ddx ) * weight_ij;
        _vy += ( dy_i + 0.5 * ddy ) * weight_ij;
      }
      else
      {
        _vx += ( dx_i + par * px_i ) * weight_ij;
        _vy += ( dy_i + par * py_i ) * weight_ij;
      }
    }
  }
 
  int rtn_flag = 0;
  if ( S_weight == 0. )
  {
    _vx = -99999.;
    _vy = -99999.;
    _vz = -99999.;
  }
  else
  {
    _vx /= S_weight;
    _vy /= S_weight;
    _vz      = -99999.;
    rtn_flag = 1;
  }
  return rtn_flag;
}
 
int FTFinder::VertexFit( int z_flag ) {
  _vx = 0.;
  _vy = 0.;
  _vz = 0.;
 
  if ( !z_flag ) { return VertexFit2D(); }
 
  if ( _tracks.size() < 2 )
  {
    _vx = -99999.;
    _vy = -99999.;
    _vz = -99999.;
    return 0;
  }
 
  FTList<float> px;
  FTList<float> py;
  FTList<float> pz;
  FTList<float> dx;
  FTList<float> dy;
  FTList<float> dz;
  FTList<float> pmag2;
  FTList<float> weight;
  FTList<float> weight_z;
  FTList<float> sigma2_r;
  FTList<float> sigma2_z;
 
  for ( auto t : _tracks )
  {
    const Lpav& la = t->lpav();
    if ( la.nc() <= 3 ) continue;
    const zav& za = t->Zav();
    if ( za.nc() > 2 && za.b() > -900 )
    {
      pmag2.push_back( 1 + za.b() * za.b() );
      pz.push_back( za.a() );
      dz.push_back( za.b() );
      sigma2_z.push_back( za.chisq() / ( za.nc() - 2 ) );
      weight_z.push_back( exp( -( za.chisq() / ( za.nc() - 2 ) ) ) );
    }
    else continue;
 
    HepVector a = la.Hpar( pivot );
 
    const float dr_i = a( 1 );
    const float px_i = -std::sin( a( 2 ) );
    const float py_i = std::cos( a( 2 ) );
 
    px.push_back( px_i );
    py.push_back( py_i );
    dx.push_back( dr_i * py_i );
    dy.push_back( -dr_i * px_i );
    sigma2_r.push_back( std::sqrt( la.chisq() ) / ( la.nc() - 3 ) );
    weight.push_back( exp( -( std::sqrt( la.chisq() ) / ( la.nc() - 3 ) ) ) );
  }
 
  if ( dx.size() < 2 )
  {
    _vx = -9999.;
    _vy = -9999.;
    _vz = -9999.;
    return 0;
  }
 
  FTList<float> vtx_x;
  FTList<float> vtx_y;
  FTList<float> vtx_z;
  FTList<float> weight2;
  FTList<float> weight2_z;
  FTList<float> vtx_chi2;
  unsigned      n_vtx( 0 );
 
  for ( int i = dx.size() - 1; i; i-- )
  {
    for ( int j = 0; j < i; j++ )
    {
      // min. chi2 fit w/ line approximation
      const float pij_x = py[i] * pz[j] - pz[i] * py[j];
      const float pij_y = pz[i] * px[j] - px[i] * pz[j];
      const float pij_z = px[i] * py[j] - py[i] * px[j];
 
      if ( pij_x == 0.0f && pij_y == 0.0f && pij_z == 0.0f ) continue;
 
      const float sr = sigma2_r[i] + sigma2_r[j];
      const float sz = sigma2_z[i] + sigma2_z[j];
 
      const float ddx = dx[i] - dx[j];
      const float ddy = dy[i] - dy[j];
      const float ddz = dz[i] - dz[j];
 
      const float pij_mag2 = pij_x * pij_x / ( sr * sz ) + pij_y * pij_y / ( sr * sz ) +
                             pij_z * pij_z / ( sr * sr );
 
      const float d_i = ( pij_x * ( py[j] * ddz - pz[j] * ddy ) / ( sr * sz ) +
                          pij_y * ( pz[j] * ddx - px[j] * ddz ) / ( sr * sz ) +
                          pij_z * ( px[j] * ddy - py[j] * ddx ) / ( sr * sr ) ) /
                        pij_mag2;
 
      const float d_j = ( pij_x * ( py[i] * ddz - pz[i] * ddy ) / ( sr * sz ) +
                          pij_y * ( pz[i] * ddx - px[i] * ddz ) / ( sr * sz ) +
                          pij_z * ( px[i] * ddy - py[i] * ddx ) / ( sr * sr ) ) /
                        pij_mag2;
 
      const float vtx_x_i = dx[i] + px[i] * d_i;
      const float vtx_y_i = dy[i] + py[i] * d_i;
      const float vtx_z_i = dz[i] + pz[i] * d_i;
 
      const float vtx_x_j = dx[j] + px[j] * d_j;
      const float vtx_y_j = dy[j] + py[j] * d_j;
      const float vtx_z_j = dz[j] + pz[j] * d_j;
 
      const float weight_ij = weight[i] + weight[j];
      vtx_x.push_back( ( weight[j] * vtx_x_i + weight[i] * vtx_x_j ) / weight_ij );
      vtx_y.push_back( ( weight[j] * vtx_y_i + weight[i] * vtx_y_j ) / weight_ij );
      vtx_z.push_back( ( weight_z[j] * vtx_z_i + weight_z[i] * vtx_z_j ) /
                       ( weight_z[i] + weight_z[j] ) );
 
      weight2.push_back( exp( -sr ) );
      weight2_z.push_back( exp( -sz ) );
      vtx_chi2.push_back( ( ( vtx_x_i - vtx_x_j ) * ( vtx_x_i - vtx_x_j ) +
                            ( vtx_y_i - vtx_y_j ) * ( vtx_y_i - vtx_y_j ) ) /
                              sr +
                          ( vtx_z_i - vtx_z_j ) * ( vtx_z_i - vtx_z_j ) / sz );
      n_vtx++;
    }
  }
 
  float S_weight( 0.0f );
  float S_weight_z( 0.0f );
  for ( int i = 0; i < vtx_chi2.size(); i++ )
  {
    if ( vtx_chi2[i] > 10. ) continue;
    float w( std::exp( -vtx_chi2[i] ) );
    _vx += vtx_x[i] * weight2[i] * w;
    _vy += vtx_y[i] * weight2[i] * w;
    _vz += vtx_z[i] * weight2_z[i] * w;
    S_weight += weight2[i] * w;
    S_weight_z += weight2_z[i] * w;
  }
 
  int rtn_flag = 0;
  if ( S_weight <= 0. )
  {
    _vx = -9999.;
    _vy = -9999.;
  }
  else
  {
    _vx /= S_weight;
    _vy /= S_weight;
    rtn_flag = 1;
  }
 
  if ( !z_flag ) return rtn_flag;
 
  if ( S_weight_z <= 0. ) { _vz = -9999.; }
  else
  {
    _vz /= S_weight_z;
    rtn_flag++;
  }
  return rtn_flag;
}
 
int FTFinder::findBestVertex() {
  if ( _tracks.size() < 2 )
  {
    _vx = -9999.;
    _vy = -9999.;
    _vz = -9999.;
    return 0;
  }
 
  float min_dr = 9999.;
  float phi0   = 9999.;
  for ( auto t : _tracks )
  {
    HepVector a = t->lpav().Hpar( pivot );
    if ( fabs( a( 1 ) ) < fabs( min_dr ) )
    {
      min_dr = a( 1 );
      phi0   = a( 2 );
    }
  }
  _vx = min_dr * cos( phi0 );
  _vy = min_dr * sin( phi0 );
  return 1;
}
 
int FTFinder::CorrectEvtTiming() {
  float weight_sum = 0.;
  float dt_sum2    = 0.;
  for ( auto t : _tracks )
  {
    float dt_sum  = 0.;
    float dtt_sum = 0.;
    int   nHits   = 0;
 
    const Lpav& la = t->lpav();
 
    const FTList<FTSegment*>& axial_sgmnts = t->axial_segments();
 
    for ( auto axial_seg : axial_sgmnts )
    {
      FTList<FTWire*>& hits = axial_seg->wireHits();
 
      for ( auto h : hits )
      {
        const float x  = h->x();
        const float y  = h->y();
        double      d0 = fabs( la.d( (double)x, (double)y ) );
        if ( d0 >= 0.47f * h->layer()->csize() ) continue;
        nHits++;
 
        float dt = x2t( h->layer(), d0 ) - h->time();
        dt_sum += dt;
        dtt_sum += ( dt * dt );
      }
    }
 
    if ( !nHits ) continue;
 
    float weight_t = exp( -( dtt_sum - ( dt_sum * dt_sum / nHits ) ) / ( nHits * 1600 ) );
    weight_sum += ( nHits * weight_t );
    dt_sum2 += ( dt_sum * weight_t );
  }
  return int( dt_sum2 / weight_sum );
}
 
#ifndef OnlineMode
void FTFinder::makeMdst( void ) {
  int Ntable( 0 );
 
  for ( auto t : _tracks )
  {
 
    FTTrack& trk = *t;
 
    const FTList<FTSegment*>& axialSegments  = trk.axial_segments();
    const FTList<FTSegment*>& stereoSegments = trk.stereo_segments();
 
    for ( auto seg : axialSegments )
    {
      FTList<FTWire*>& wires = seg->wireHits();
      for ( auto w : wires ) { g_ncell->fill( w->getWireId() ); }
    }
 
    for ( auto seg : stereoSegments )
    {
      FTList<FTWire*>& wires = seg->wireHits();
      for ( auto w : wires ) { g_ncell->fill( w->getWireId() ); }
    }
 
    const Helix& h = *trk.helix();
    if ( h.tanl() < -9000. ) continue;
 
    Ntable++;
    Hep3Vector p( h.momentum() );
    HepPoint3D v( h.x() );
    // float m_charge = (h.kappa() > 0) ? 1. : -1.;
    float m_px        = p.x();
    g_px[Ntable - 1]  = p.x();
    float m_py        = p.y();
    g_py[Ntable - 1]  = p.y();
    float m_pz        = p.z();
    g_pz[Ntable - 1]  = p.z();
    g_p[Ntable - 1]   = p.mag();
    g_phi[Ntable - 1] = atan2( m_py, m_px );
 
    g_dr[Ntable - 1]    = h.dr();
    g_phi0[Ntable - 1]  = h.phi0();
    g_kappa[Ntable - 1] = h.kappa();
    g_pt[Ntable - 1]    = 1 / fabs( h.kappa() );
    g_dz[Ntable - 1]    = h.dz();
    g_tanl[Ntable - 1]  = h.tanl();
    g_theta[Ntable - 1] = atan2( 1 / fabs( h.kappa() ), (double)( m_pz ) );
    g_vx[Ntable - 1]    = v( 0 );
    g_vy[Ntable - 1]    = v( 1 );
    g_vz[Ntable - 1]    = v( 2 );
  }
  g_ntrk = Ntable;
}
#endif
 
//-------------------------------
// output tracking results to TDS by wangdy
//-------------------------------
StatusCode FTFinder::makeTds() {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
  MsgStream log( msgSvc, "FTFinder" );
 
  IDataProviderSvc* eventSvc;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc ).ignore();
 
  if ( !eventSvc )
  {
    log << MSG::FATAL << "makeTds:Could not find eventSvc" << endmsg;
    return StatusCode::FAILURE;
  }
 
  RecMdcTrackCol* trkcol = new RecMdcTrackCol;
  RecMdcHitCol*   hitcol = new RecMdcHitCol;
 
  // make RecMdcTrackCol
#ifndef OnlineMode
  log << MSG::DEBUG << "beginning to make RecMdcTrackCol" << endmsg;
#endif
 
  int trkid = 0, i = -1;
  for ( auto fttrk : tracks() )
  {
    i++;
 
    RecMdcTrack* trk = new RecMdcTrack;
 
    if ( fttrk->helix()->tanl() < -9000. )
    {
      log << MSG::DEBUG << "deleted trackId = " << i + 1
          << " due to tanl = " << fttrk->helix()->tanl() << endmsg;
      delete trk;
      continue;
    }
 
    HepPoint3D pos    = fttrk->helix()->pivot();
    int        charge = -1;
    HepVector  m_a( 5, 0 );
    m_a                   = fttrk->helix()->a();
    m_a[2]                = -m_a[2];
    HepSymMatrix m_ea     = fttrk->helix()->Ea();
    float        fiterm   = fttrk->lpav().phi( 77.0 );
    float        chi2lpav = fttrk->lpav().chisq();
    float        chi2zav  = fttrk->Zav().chisq();
 
    m_ea[0][0] = 0.0085;
    m_ea[1][1] = 0.000011;
    m_ea[2][2] = 0.0018;
    m_ea[3][3] = 1.2;
    m_ea[4][4] = 0.00026;
    m_ea[1][0] = m_ea[0][1] = -0.00029;
    m_ea[2][0] = m_ea[0][2] = charge * 0.004;
    m_ea[2][1] = m_ea[1][2] = charge * 0.00012;
    m_ea[3][0] = m_ea[0][3] = -0.017;
    m_ea[3][1] = m_ea[1][3] = 0.0;
    m_ea[3][2] = m_ea[2][3] = 0.0;
    m_ea[4][0] = m_ea[0][4] = 0.0;
    m_ea[4][1] = m_ea[1][4] = 0.0;
    m_ea[4][2] = m_ea[2][4] = 0.0;
    m_ea[4][3] = m_ea[3][4] = -0.018;
 
    trk->setTrackId( trkid );
    trk->setHelix( m_a );
    trk->setPivot( pos );
    trk->setError( m_ea );
    trk->setFiTerm( fiterm );
    trk->setChi2( chi2lpav + chi2zav );
 
#ifndef OnlineMode
    log << MSG::DEBUG << " trackId = " << i << endmsg;
    log << MSG::DEBUG << "fast-tracking  kappa " << m_a[2] << "   fast-tracking tanl "
        << m_a[4] << endmsg;
    log << MSG::DEBUG << "push_backed kappa " << trk->helix( 2 ) << "   push_backed tanl "
        << trk->helix( 4 ) << endmsg;
 
    log << MSG::DEBUG << "beginning to make hitlist and RecMdcHitCol " << endmsg;
#endif
 
    HitRefVec hitrefvec;
 
    int hitindex = 0;
 
    // axial segments part
    const FTList<FTSegment*>& seglist_ax = fttrk->axial_segments();
 
    int ntrackhits = 0;
    for ( auto seg_ax_i : seglist_ax )
    {
      FTList<FTWire*>& hitlist_ax = seg_ax_i->wireHits();
      ntrackhits += hitlist_ax.size();
 
      for ( auto wire_ax : hitlist_ax )
      {
        double           dd   = wire_ax->distance();
        double           adc  = RawDataUtil::MdcChargeChannel( wire_ax->getAdc() );
        double           tdc  = wire_ax->time();
        double           chi2 = wire_ax->getChi2();
        const Identifier mdcid =
            MdcID::wire_id( wire_ax->layer()->layerId(), wire_ax->localId() );
 
        RecMdcHit* hit = new RecMdcHit;
        hit->setId( hitindex );
        hit->setDriftDistLeft( dd );
        hit->setDriftDistRight( dd );
        hit->setAdc( adc );
        hit->setTdc( tdc );
        hit->setMdcId( mdcid );
        hit->setChisqAdd( chi2 );
 
#ifndef OnlineMode
        log << MSG::DEBUG << "Hit DriftDistLeft " << hit->getDriftDistLeft() << endmsg;
        log << MSG::DEBUG << "MDC Hit ADC " << hit->getAdc() << endmsg;
        log << MSG::DEBUG << "Hit MDC  Identifier " << hit->getMdcId() << endmsg;
        log << MSG::DEBUG << "Hit Chisq " << hit->getChisqAdd() << endmsg;
#endif
 
        hitcol->push_back( hit );
        SmartRef<RecMdcHit> refhit( hit );
        hitrefvec.push_back( refhit );
      }
    }
 
    // stereo segments part
    const FTList<FTSegment*>& seglist_st = fttrk->stereo_segments();
 
    int ntrackster = 0;
 
    for ( auto seg_st_i : seglist_st )
    {
      FTList<FTWire*>& hitlist_st = seg_st_i->wireHits();
 
      ntrackhits += hitlist_st.size();
      ntrackster += hitlist_st.size();
 
      for ( auto wire_st : hitlist_st )
      {
        double           dd   = wire_st->distance();
        double           adc  = RawDataUtil::MdcChargeChannel( wire_st->getAdc() );
        double           tdc  = wire_st->time();
        double           chi2 = wire_st->getChi2();
        const Identifier mdcid =
            MdcID::wire_id( wire_st->layer()->layerId(), wire_st->localId() );
 
        RecMdcHit* hit = new RecMdcHit;
        hit->setId( hitindex );
        hit->setDriftDistLeft( dd );
        hit->setDriftDistRight( dd );
        hit->setAdc( adc );
        hit->setTdc( tdc );
        hit->setMdcId( mdcid );
        hit->setChisqAdd( chi2 );
 
#ifndef OnlineMode
        log << MSG::DEBUG << "Z Hit DriftDistLeft " << hit->getDriftDistLeft() << endmsg;
        log << MSG::DEBUG << "Z MDC Hit ADC " << hit->getAdc() << endmsg;
        log << MSG::DEBUG << "Z Hit MDC  Identifier " << hit->getMdcId() << endmsg;
        log << MSG::DEBUG << "Z Hit Chisq " << hit->getChisqAdd() << endmsg;
#endif
 
        hitcol->push_back( hit );
        SmartRef<RecMdcHit> refhit( hit );
        hitrefvec.push_back( refhit );
      }
    }
 
    trk->setNhits( ntrackhits );
    trk->setNster( ntrackster );
    trk->setVecHits( hitrefvec );
    trkcol->push_back( trk );
    trkid++;
 
#ifndef OnlineMode
    g_naxialhit->fill( ( ntrackhits - ntrackster ), 1.0 );
    g_nstereohit->fill( ntrackster, 1.0 );
    g_nhit->fill( ntrackhits, 1.0 );
#endif
  }
 
  IDataManagerSvc* dataManSvc = dynamic_cast<IDataManagerSvc*>( eventSvc );
 
  DataObject* aTrackCol;
  eventSvc->findObject( "/Event/Recon/RecMdcTrackCol", aTrackCol ).ignore();
  if ( aTrackCol )
  {
    dataManSvc->clearSubTree( "/Event/Recon/RecMdcTrackCol" ).ignore();
    eventSvc->unregisterObject( "/Event/Recon/RecMdcTrackCol" ).ignore();
  }
 
  DataObject* aRecHitCol;
  eventSvc->findObject( "/Event/Recon/RecMdcHitCol", aRecHitCol ).ignore();
  if ( aRecHitCol )
  {
    dataManSvc->clearSubTree( "/Event/Recon/RecMdcHitCol" ).ignore();
    eventSvc->unregisterObject( "/Event/Recon/RecMdcHitCol" ).ignore();
  }
 
  // register RecMdcHitCol
  StatusCode hitsc;
  hitsc = eventSvc->registerObject( "/Event/Recon/RecMdcHitCol", hitcol );
  if ( !hitsc.isSuccess() )
  {
    log << MSG::FATAL << "Could not register RecMdcHitCol" << endmsg;
    return hitsc;
  }
 
#ifndef OnlineMode
  log << MSG::DEBUG << "RecMdcHitCol registered successfully!" << endmsg;
 
  RecMdcTrackCol::iterator bef = trkcol->begin();
  for ( ; bef != trkcol->end(); bef++ )
  {
    log << MSG::DEBUG << " registered kappa" << ( *bef )->helix( 2 ) << "   registered tanl"
        << ( *bef )->helix( 4 ) << endmsg;
  }
#endif
 
  // register RecMdcTrackCol
  StatusCode trksc;
  trksc = eventSvc->registerObject( "/Event/Recon/RecMdcTrackCol", trkcol );
  if ( !trksc.isSuccess() )
  {
    log << MSG::FATAL << "Could not register RecMdcTrackCol" << endmsg;
    return trksc;
  }
 
#ifndef OnlineMode
  log << MSG::DEBUG << "RecMdcTrackCol registered successfully!" << endmsg;
 
  // check the result:RecMdcHitCol
  SmartDataPtr<RecMdcHitCol> newhitCol( eventSvc, "/Event/Recon/RecMdcHitCol" );
  if ( !newhitCol )
  {
    log << MSG::FATAL << "Could not find RecMdcHit" << endmsg;
    return ( StatusCode::FAILURE );
  }
 
  log << MSG::DEBUG << "Begin to check RecMdcHitCol" << endmsg;
 
  RecMdcHitCol::iterator iter_hit = newhitCol->begin();
  for ( ; iter_hit != newhitCol->end(); iter_hit++ )
  {
    log << MSG::DEBUG << "retrieved MDC Hit:"
        << "   DDL " << ( *iter_hit )->getDriftDistLeft() << "   DDR "
        << ( *iter_hit )->getDriftDistRight() << "   ADC " << ( *iter_hit )->getAdc()
        << endmsg;
  }
 
  // check the result:RecMdcTrackCol
  SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc, "/Event/Recon/RecMdcTrackCol" );
  if ( !newtrkCol )
  {
    log << MSG::FATAL << "Could not find RecMdcTrackCol" << endmsg;
    return ( StatusCode::FAILURE );
  }
 
  log << MSG::DEBUG << "Begin to check RecMdcTrackCol" << endmsg;
  RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
  for ( ; iter_trk != newtrkCol->end(); iter_trk++ )
  {
    log << MSG::DEBUG << "retrieved MDC track:"
        << "Track Id: " << ( *iter_trk )->trackId() << "   Pivot: " << ( *iter_trk )->poca()[0]
        << " " << ( *iter_trk )->poca()[1] << " " << ( *iter_trk )->poca()[2] << endmsg;
 
    log << MSG::DEBUG << "hitList of this  track:" << endmsg;
    HitRefVec           gothits   = ( *iter_trk )->getVecHits();
    HitRefVec::iterator it_gothit = gothits.begin();
    for ( ; it_gothit != gothits.end(); it_gothit++ )
    {
      log << MSG::DEBUG << "hits Id: " << ( *it_gothit )->getId() << "   hits DDL&DDR "
          << ( *it_gothit )->getDriftDistLeft() << "   hits MDC IDentifier "
          << ( *it_gothit )->getMdcId() << endmsg << "   hits TDC " << ( *it_gothit )->getTdc()
          << "   hits ADC " << ( *it_gothit )->getAdc() << endmsg;
    }
  }
#endif
 
  return StatusCode::SUCCESS;
}
 
void FTFinder::term() {
  for ( auto t : _tracks ) delete t;
  for ( auto s : _linkedSegments ) delete s;
 
  if ( param->_doIt ) clear();
}