MdcTrackList.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: MdcTrackList.cxx,v 1.25 2017/02/22 05:53:33 zhangy Exp $
//
// Environment:
//      Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author(s):
//  Steve Schaffner
//  Zhang Xueyao(zhangxy@hepg.sdu.edu.cn) Migration for BESIII
//  Zhang Yao(zhangyao@ihep.ac.cn)
//
//--------------------------------------------------------------------------
#include "MdcTrkRecon/MdcTrackList.h"
#include "CLHEP/Alist/AList.h"
#include "MdcData/MdcHit.h"
#include "MdcData/MdcHitMap.h"
#include "MdcData/MdcRecoHitOnTrack.h"
#include "MdcGeom/BesAngle.h"
#include "MdcGeom/Constants.h"
#include "MdcGeom/MdcDetector.h"
#include "MdcRawEvent/MdcDigi.h"
#include "MdcTrkRecon/GmsList.h"
#include "MdcTrkRecon/MdcSeg.h"
#include "MdcTrkRecon/MdcSegGrouperAx.h"
#include "MdcTrkRecon/MdcSegGrouperSt.h"
#include "MdcTrkRecon/MdcSegInfo.h"
#include "MdcTrkRecon/MdcSegInfoAxialO.h"
#include "MdcTrkRecon/MdcSegList.h"
#include "MdcTrkRecon/MdcTrack.h"
#include "MdcTrkRecon/MdcTrackParams.h"
#include "MdcTrkRecon/mdcWrapWire.h"
#include "PatBField/BField.h"
#include "TrkBase/TrkErrCode.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrkBase/TrkFit.h"
#include "TrkBase/TrkFitStatus.h"
#include "TrkBase/TrkHitOnTrk.h"
#include "TrkBase/TrkHotList.h"
#include "TrkBase/TrkRecoTrk.h"
#include "TrkFitter/TrkCircleMaker.h"
#include "TrkFitter/TrkHelixMaker.h"
#include <iostream>
#include <math.h>
 
// yzhang hist cut
#include "AIDA/IHistogram1D.h"
#include "GaudiKernel/NTuple.h"
#include "Rtypes.h"
extern int                 haveDigiTk[43][288];
extern double              haveDigiDrift[43][288];
extern AIDA::IHistogram1D* g_cirTkChi2;
extern AIDA::IHistogram1D* g_3dTkChi2;
extern AIDA::IHistogram1D* g_fitNAct;
extern AIDA::IHistogram1D* g_pickHitWire;
// extern AIDA::IHistogram2D*  g_predDriftVsMcTk;
extern NTuple::Tuple*        m_tuplePick;
extern NTuple::Item<long>    m_pickIs2D;
extern NTuple::Item<long>    m_pickLayer;
extern NTuple::Item<long>    m_pickNCell;
extern NTuple::Item<long>    m_pickNCellWithDigi;
extern NTuple::Array<long>   m_pickWire;
extern NTuple::Array<double> m_pickPredDrift;
extern NTuple::Array<double> m_pickDrift;
extern NTuple::Array<double> m_pickDriftTruth;
extern NTuple::Array<int>    m_pickPhizOk;
extern NTuple::Array<double> m_pickZ;
extern NTuple::Array<double> m_pickResid;
extern NTuple::Array<double> m_pickSigma;
extern NTuple::Array<double> m_pickPhiLowCut;
extern NTuple::Array<double> m_pickPhiHighCut;
extern NTuple::Array<double> m_pickDriftCut;
extern NTuple::Array<long>   m_pickMcTk;
extern NTuple::Array<double> m_pickPt;
extern NTuple::Array<double> m_pickCurv;
 
// zhangy hist cut
 
#ifdef MDCPATREC_TIMETEST
// tau include
#  include <Profile/Profiler.h>
#  include <pthread.h>
#  include <sys/time.h>
#endif
// zhangy hist
using MdcPatRec::BField;
using std::cout;
using std::endl;
 
// *************************************************************************
MdcTrackList::MdcTrackList( const MdcTrackParams& tkPar ) : MdcTrackListBase( tkPar ) {
  // *************************************************************************
  return;
}
 
//------------------------------------------------------------------------
MdcTrackList::~MdcTrackList() {}
//------------------------------------------------------------------------
 
// *************************************************************************
int MdcTrackList::createFromSegs( MdcSegList* segs, const MdcHitMap* map,
                                  const MdcDetector* gm, TrkContext& context,
                                  double bunchTime ) {
  // *************************************************************************
  // Forms tracks out of list of superlayer segments
  int nTracks = 0;
 
  m_debug = tkParam.lPrint; // yzhang debug
 
  // Create empty list of stereo segs (to save time)
  MdcSegGrouperSt stereoSegs( gm, tkParam.lPrint );
  stereoSegs.setNoInner( m_noInner ); // yzhang add 2016-10-12
 
  // *** Create r-phi track
 
  // Create list of axial segments, treated as on tracks from origin
#ifdef MDCPATREC_TIMETEST
  TAU_PROFILE_TIMER( timer8, "fill ax seg", "int ()", TAU_DEFAULT );
  TAU_PROFILE_START( timer8 );
#endif
  MdcSegGrouperAx origSegs( gm, tkParam.lPrint );
  origSegs.fillWithSegs( segs );
  // std::cout << "Plot segs after ax fillWithSegs " << std::endl;//yzhang debug
  // segs->plot(0);//yzhang debug
#ifdef MDCPATREC_TIMETEST
  TAU_PROFILE_STOP( timer8 );
#endif
  MdcSeg* seed;
  bool    goodOnly = true;
  bool    useBad   = ( segs->count() < 400 ); // if true, use non-good seeds
  // bool useBad = false;
  segs->resetSeed( gm );
  MdcTrack* trialTrack = 0;
 
  while ( 1 )
  {
    seed = segs->getSeed( 0, goodOnly );
    if ( seed == 0 && goodOnly && useBad )
    {
      segs->resetSeed( gm );
      goodOnly = false;
      continue;
    }
    else if ( seed == 0 && ( !goodOnly || !useBad ) ) { break; }
 
    if ( 3 == m_debug ) dumpSeed( seed, goodOnly ); // yzhang debug
 
    // Pass through first superlayer of MDC required, Field layer No5 = 12.135 cm
    if ( fabs( ( (MdcSegInfoAxialO*)seed->info() )->curv() ) > 0.0417 )
      continue; // curv cut yzhang
    delete trialTrack;
    trialTrack = 0;
 
    //---------Combine Ax segs--------
#ifdef MDCPATREC_TIMETEST
    TAU_PROFILE_TIMER( timer3, "combine ax seg", "int ()", TAU_DEFAULT );
    TAU_PROFILE_START( timer3 );
#endif
    int success =
        origSegs.combineSegs( trialTrack, seed, context, bunchTime, tkParam.maxSegChisqO );
#ifdef MDCPATREC_TIMETEST
    TAU_PROFILE_STOP( timer3 );
#endif
 
    if ( 3 <= m_debug )
    {
      cout << " ***** Ax.combineSegs success?  " << success << "\n";
      dumpAxCombine( trialTrack ); // yzhang debug
    }
    if ( !success ) continue;
 
      //--------Finish circle--------
#ifdef MDCPATREC_TIMETEST
    TAU_PROFILE_TIMER( timer4, "circle fitting", "int ()", TAU_DEFAULT );
    TAU_PROFILE_START( timer4 );
#endif
    success = finishCircle( *trialTrack, map, gm );
#ifdef MDCPATREC_TIMETEST
    TAU_PROFILE_STOP( timer4 );
#endif
 
    if ( 3 <= m_debug )
    {
      cout << "finishCircle success?  " << success << "\n";
      dumpCircle( trialTrack ); // yzhang debug
    }
 
    if ( !success ) continue;
    //--------Make sure it really did come from origin--------
    const TrkFit* tkFit = trialTrack->track().fitResult();
    assert( tkFit != 0 );
    TrkExchangePar par = tkFit->helix( 0.0 );
    // dumpD0(par);
 
    //------------------d0 cut-------------------
    // 2010-03-31 , m_d0Cut from 0 to epsilon
 
    // std::cout<<  __FILE__ <<"  d0 cut------------"<< fabs(par.d0())<<" d0cut "<< m_d0Cut <<
    // std::endl;
    if ( ( m_d0Cut > Constants::epsilon ) && ( fabs( par.d0() ) > m_d0Cut ) )
    {
      if ( tkParam.lPrint > 1 )
        cout << __FILE__ << " Killing track by d0:" << par.d0() << ">" << m_d0Cut << endl;
      continue;
    }
 
    //------------------pt cut-------------------
    if ( ( fabs( m_ptCut ) > 0. ) && ( fabs( tkFit->pt( 0. ) ) < m_ptCut ) )
    {
      if ( tkParam.lPrint > 1 )
        cout << __FILE__ << " Killing track by pt:" << tkFit->pt( 0. ) << "<" << m_ptCut
             << endl;
      continue;
    }
 
    // *** End r-phi track-finding
    if ( 3 <= m_debug ) std::cout << " *** End r-phi track-finding " << std::endl;
 
      //--------Add stereo to taste--------
#ifdef MDCPATREC_TIMETEST
    TAU_PROFILE_TIMER( timer5, "fill st seg", "int ()", TAU_DEFAULT );
    TAU_PROFILE_START( timer5 );
#endif
    stereoSegs.fillWithSegs( segs, trialTrack );
#ifdef MDCPATREC_TIMETEST
    TAU_PROFILE_STOP( timer5 );
#endif
    if ( 3 <= m_debug )
    {
      // dumpStFill();//yzhang debug
      std::cout << std::endl << "----------------------------------------" << std::endl;
      std::cout << " Segment list after St.fillWithSegs " << std::endl;
      stereoSegs.dumpSegList();
    }
 
#ifdef MDCPATREC_TIMETEST
    TAU_PROFILE_TIMER( timer6, "combine st seg", "int ()", TAU_DEFAULT );
    TAU_PROFILE_START( timer6 );
#endif
    success = stereoSegs.combineSegs( trialTrack, 0, context, bunchTime, tkParam.maxSegChisqO,
                                      tkParam.combineByFitHits );
#ifdef MDCPATREC_TIMETEST
    TAU_PROFILE_STOP( timer6 );
#endif
 
    if ( 3 <= m_debug )
    {
      cout << " ***** St.combineSegs success?  " << success << "\n";
      dumpStCombine( trialTrack );
    }
 
    //--------Finish 3-d track--------
    if ( success )
    {
#ifdef MDCPATREC_TIMETEST
      TAU_PROFILE_TIMER( timer7, "helix fitting", "int ()", TAU_DEFAULT );
      TAU_PROFILE_START( timer7 );
#endif
      success = finishHelix( *trialTrack, map, gm );
#ifdef MDCPATREC_TIMETEST
      TAU_PROFILE_STOP( timer7 );
#endif
    } // end if (success -- st)
 
    if ( 3 == m_debug ) { dumpHelix( trialTrack ); }
    if ( !success ) continue;
 
    //------------------d0 cut after helix fitting-------------------
    // yzhang add 2011-08-01
    double d0par = trialTrack->track().fitResult()->helix( 0. ).d0();
    if ( ( m_d0Cut > Constants::epsilon ) && ( fabs( d0par ) > m_d0Cut ) )
    {
      if ( tkParam.lPrint > 1 )
      {
        cout << __FILE__ << " Killing track by d0 after 3d fit:" << d0par << ">" << m_d0Cut
             << endl;
      }
      continue;
    }
 
    double z0par = trialTrack->track().fitResult()->helix( 0. ).z0();
    if ( ( m_z0Cut > Constants::epsilon ) && ( fabs( z0par ) > m_z0Cut ) )
    {
      if ( tkParam.lPrint > 1 )
      { cout << __FILE__ << " Killing track by z0:" << z0par << ">" << m_z0Cut << endl; }
      continue;
    }
 
    nTracks++;
    append( trialTrack ); // Add to list of Tracks
 
    trialTrack = 0;
 
    if ( 3 == m_debug ) std::cout << " ***** End one track-finding *****" << std::endl;
  } // end while(1)
 
  delete trialTrack;
  return nTracks;
}
 
/***************************************************************************/
int MdcTrackList::pickHits( MdcTrack* trk, const MdcHitMap* map, const MdcDetector* gm,
                            bool pickAm ) {
 
  /***************************************************************************/
 
  // Selects candidate hits along track;
  // sorts them into "active" (small residual) and "inactive" (large resid);
  // throws away hits separated from main group by large gaps. //?? FIXME
  // pickAm => pick the ambiguity for hits already on track
  // Return # of active hits found
 
  bool is2d = trk->track().status()->is2d();
  if ( 6 == tkParam.lPrint )
  {
    cout << "Before pickHits";
    if ( is2d ) cout << " 2d:";
    else cout << " 3d:";
    cout << endl;
  }
 
  int                     nFound = 0;
  int                     nCand  = 0;        // cells tried
  double                  rMin, rMax;        // min & max search radii for this track
  double                  rEnter, rExit;     // radii at which track enters, exits layer
  BesAngle                phiEnter, phiExit; // yzhang change
  int                     wireLow, wireHigh;
  double                  phiLow, phiHigh;
  int                     nCell;
  MdcHit*                 thisHit;
  HepAList<MdcHitOnTrack> localHistory; // temporary list of added hits
  // until bogus hits are chucked
  double cellWidth;
  double goodDriftCut; // missing hits with predDrift > goodDriftCut don't count in figuring
                       // gaps
  double aresid        = 0.; // = abs(resid)
  int    firstInputHit = -1; // first hit in firstInputLayer/lastInputLayer region
  int    lCurl         = 0;  // reached curl point
 
  //****************************************************/
  const MdcLayer* firstInputLayer = trk->firstLayer();
  const MdcLayer* lastInputLayer  = trk->lastLayer();
 
  double        bunchTime = trk->track().trackT0();
  const TrkFit* tkFit     = trk->track().fitResult();
  assert( tkFit != 0 );
  const TrkFitStatus* tkStatus = trk->track().status();
  assert( tkStatus != 0 );
  TrkHitList* hitList = trk->track().hits();
  assert( hitList != 0 );
  TrkExchangePar par     = tkFit->helix( 0.0 );
  double         d0      = par.d0();
  double         curv    = 0.5 * par.omega(); //!!! change to using omega itself
  double         sinPhi0 = sin( par.phi0() );
  double         cosPhi0 = cos( par.phi0() );
 
  // *** Set min and max radius for track
  rMin         = gm->firstLayer()->rIn();
  double absd0 = fabs( d0 );
  if ( absd0 > rMin ) rMin = absd0 + Constants::epsilon;
 
  bool   willCurl = false;
  double rCurl    = fabs( d0 + 1. / curv );
  rMax            = gm->lastLayer()->rOut();
  // std::cout<<  __FILE__ <<"  "<< __LINE__<<" rCurl "<< rCurl <<" rMax "<< rMax << std::endl;
  if ( rCurl < rMax )
  {
    willCurl = true;
    rMax     = rCurl - Constants::epsilon;
  }
  // std::cout<<" willCurl "<< willCurl  << std::endl;
  bool reachedLastInput = false;
  bool hasCurled        = false; // hit found past curl point
 
  // yzhang add skip layer with hot, 2011-05-03
  bool isHotOnLayer[43];
  if ( tkParam.pickSkipExistLayer )
  {
    for ( int ii = 0; ii < 43; ii++ ) isHotOnLayer[ii] = false;
    for ( TrkHitList::hot_iterator ihit( hitList->begin() ); ihit != hitList->end(); ++ihit )
    { isHotOnLayer[ihit->layerNumber()] = true; }
  }
 
  // *** Loop through layers in view, looking for the hits
  // assumes axial inner
  for ( const MdcLayer* layer = gm->firstLayer(); layer != 0;
        layer =
            ( ( !lCurl ) ? ( ( hasCurled ) ? gm->prevLayer( layer ) : gm->nextLayer( layer ) )
                         : layer ) )
  {
 
    if ( lCurl )
    {
      lCurl     = 0;
      hasCurled = true;
    }
    if ( tkStatus->is2d() && layer->view() != 0 ) continue;
 
    if ( tkParam.pickSkipExistLayer && isHotOnLayer[layer->layNum()] )
      continue; // yzhang add 2011-05-03
 
    // std::cout<<  __FILE__ <<"  "<< __LINE__  << " lCurl "<< lCurl
    //<<" hasCurled "<< hasCurled <<" layer "<< layer->layNum() << std::endl;
    bool lContinue = true;
    if ( tkParam.pickUitlLastLayer )
    { // yzhang change 2010-09-10
      if ( layer == lastInputLayer ) reachedLastInput = true;
    }
 
    // *** Find enter and exit points
    if ( hasCurled )
    {
      rEnter = layer->rOut();
      if ( rEnter < rMin )
      {
        // std::cout<<  __FILE__ <<"  "<< __LINE__
        //<<" rEnter<rMin "<<rEnter<<" "<<rMin<< std::endl;
        break;
      }
      rExit = layer->rIn();
      if ( rExit < rMin ) rExit = rMin;
      if ( rEnter > rMax ) rEnter = rMax;
 
      // std::cout<<  __FILE__ <<" hasCurled: rEnter  "<< rEnter
      //<<" rExit " << rExit << " rMin "<<rMin<<" rMax "<<rMax<< std::endl;
    }
    else
    {
      rEnter = layer->rIn();
      rExit  = layer->rOut();
      // std::cout<<  __FILE__ <<" NOT hasCurled: rEnter  "<< rEnter
      //<<" rExit " << rExit << " rMin "<<rMin<<" rMax "<<rMax<< std::endl;
      if ( rExit < rMin ) continue;
 
      if ( willCurl )
      {
        if ( rEnter > rMax )
        {
          // std::cout<<  __FILE__<<" Reached curl point before entering layer"<< std::endl;
          //  Reached curl point before entering layer
          hasCurled = 1;
          continue;
        }
        if ( rExit > rMax )
        {
          lCurl = 1;
          rExit = rMax;
        }
      }
      else
      { // not a potential curler
        // std::cout<<  __FILE__ <<"  "<< __LINE__  <<" not a potential curler"<< std::endl;
        if ( rEnter > rMax )
        {
          // std::cout<<  __FILE__ <<" rEnter> rMax "<< rEnter << std::endl;
          break;
        }
        if ( rExit > rMax ) rExit = rMax;
      }
    } // end if curled already
 
    nCell     = layer->nWires();
    cellWidth = Constants::twoPi * layer->rMid() / nCell; //??
    // don't count outer xmm of cell
    goodDriftCut = 0.5 * cellWidth * M_SQRT2 + tkParam.pickHitMargin;
    // goodDriftCut = 0.5 * cellWidth - tkParam.pickHitMargin;//yzhang change 2012-08-17
    double deltaPhiCellWidth = 0.5 * ( cellWidth * M_SQRT2 ) / layer->rMid();
 
    //**** Find entrance and exit phi's
    BesAngle phiTemp( 0.0 );
    int      ierror = trk->projectToR( rEnter, phiTemp, hasCurled );
    phiEnter        = phiTemp.posRad();
    if ( ierror != 0 )
    {
      if ( 6 == tkParam.lPrint )
        std::cout << " ErrMsg(warning) "
                  << "Error in MdcTrackList::pickHits projection, ierror " << ierror
                  << std::endl;
      continue; // yzhang 2011-04-14
    }
    ierror  = trk->projectToR( rExit, phiTemp, hasCurled );
    phiExit = phiTemp.posRad();
    if ( ierror != 0 )
    {
      if ( 6 == tkParam.lPrint )
        std::cout << " ErrMsg(warning) "
                  << "Error in MdcTrackList::pickHits projection, ierror " << ierror
                  << std::endl;
      continue; // yzhang 2011-04-14
    }
 
    if ( 6 == tkParam.lPrint )
    {
      std::cout << endl << "--pickHit of layer " << layer->layNum() << "--" << std::endl;
      std::cout << "  track phiEnter " << phiEnter.deg() << " phiExit " << phiExit.deg()
                << " degree" << std::endl;
      std::cout << "  cell width " << 360. / layer->nWires() << " dPhiz "
                << layer->dPhiz() * Constants::radToDegrees << " deltaPhiCellWidth "
                << 0.5 * ( cellWidth * M_SQRT2 ) / layer->rMid() * Constants::radToDegrees
                << std::endl;
      std::cout << "  maxInactiveResid " << tkParam.maxInactiveResid << " pickHitPhiFactor "
                << tkParam.pickHitPhiFactor << std::endl;
      std::cout << "  goodDriftCut " << goodDriftCut << "=sqrt(2)*0.5*" << cellWidth << "+"
                << tkParam.pickHitMargin << std::endl;
    }
 
    double Bz = trk->track().bField().bFieldZ();
    // std::cout<<  " facotr "<<tkParam.pickHitPhiFactor<< " dPhiz "<<layer->dPhiz()
    //<< " factor*dPhiz "<<layer->dPhiz()*tkParam.pickHitPhiFactor<< std::endl;
    double t_phiHit = -999.;
    if ( curv * Bz <= 0.0 )
    {
      // positive track in minus Bz
      phiLow  = phiEnter;
      phiHigh = phiExit;
      // Allow some slop in phi
      phiLow -= tkParam.maxInactiveResid / rEnter;
      phiHigh += tkParam.maxInactiveResid / rExit;
    }
    else
    {
      phiLow  = phiExit;
      phiHigh = phiEnter;
      phiLow -= tkParam.maxInactiveResid / rExit;
      phiHigh += tkParam.maxInactiveResid / rEnter;
    }
    // yzhang fix cross x axis bug 2011-04-10
    if ( ( phiHigh > 0 && phiLow < 0 ) ) { phiLow += Constants::twoPi; }
    else if ( ( phiHigh < 0 && phiLow > 0 ) ) { phiHigh += Constants::twoPi; }
 
    double lowFloat  = nCell / Constants::twoPi * ( phiLow - layer->phiOffset() ) + 0.5;
    double highFloat = nCell / Constants::twoPi * ( phiHigh - layer->phiOffset() ) + 0.5;
    wireLow          = ( lowFloat >= 0.0 ) ? int( lowFloat ) : int( floor( lowFloat ) );
    wireHigh         = ( highFloat >= 0.0 ) ? int( highFloat ) : int( floor( highFloat ) );
 
    if ( 6 == tkParam.lPrint )
    {
      std::cout << "  wireLow " << wireLow << " wireHigh " << wireHigh << " phiLow "
                << phiLow * 180. / Constants::pi << " phiHigh "
                << phiHigh * 180. / Constants::pi << std::endl;
    }
    // *** Loop through the predicted hit wires
    int tempDiff = 0;
    if ( g_pickHitWire )
    {
      int tempN = Constants::maxCell[layer->layNum()];
      if ( wireLow > tempN / 2. && wireHigh < tempN / 2. )
      {
        g_pickHitWire->fill( wireHigh + tempN - wireLow );
        tempDiff = wireHigh + tempN - wireLow + 1;
      }
      else
      {
        g_pickHitWire->fill( wireHigh - wireLow );
        tempDiff = wireHigh - wireLow + 1;
      }
    } // yzhang hist cut
 
    if ( wireLow > wireHigh ) wireLow -= nCell; // yzhang 2011-05-16
    long t_iHit = 0;
    for ( int thisWire = wireLow; thisWire <= wireHigh; thisWire++ )
    {
      int corrWire = mdcWrapWire( thisWire, nCell );
      thisHit      = map->hitWire( layer->layNum(), corrWire );
 
      if ( 6 == tkParam.lPrint )
      {
        if ( thisHit != 0 )
        {
          cout << endl << "test hit ";
          thisHit->print( std::cout );
        }
        else
          std::cout << endl
                    << "test hit (" << layer->layNum() << "," << corrWire << ")" << std::endl;
      }
 
      double tof       = 0.;
      double z         = 0.;
      double driftDist = 0.;
      // Calculate expected drift distance
      double               delx, dely;
      double               resid = 0., predDrift = 0.;
      int                  ambig  = 0;
      const MdcHitOnTrack* alink  = 0;
      double               mcTkId = -9999; // yzhang for tuple 2011-06-28
      if ( thisHit == 0 )
      {
        delx      = -d0 * sinPhi0 - layer->xWire( corrWire );
        dely      = d0 * cosPhi0 - layer->yWire( corrWire );
        predDrift = cosPhi0 * dely - sinPhi0 * delx + curv * ( delx * delx + dely * dely );
        if ( 6 == tkParam.lPrint ) cout << "No hit. predDrift=" << predDrift << endl;
        continue;
      }
      else
      {
        if ( m_tuplePick ) mcTkId = thisHit->digi()->getTrackIndex();
        // look for existing hit
        if ( thisHit->getHitOnTrack( &( trk->track() ) ) == 0 )
        {
          alink = 0; // yzhang temp
        }
        else
        {
          alink = thisHit->getHitOnTrack( &( trk->track() ) )->mdcHitOnTrack(); // yzhang temp
          if ( 6 == tkParam.lPrint ) { cout << " existing hot; "; }
        }
 
        if ( alink == 0 || pickAm )
        {
          if ( ( !tkStatus->is2d() ) && layer->view() != 0 )
          {
            double rc = 9999.;
            if ( fabs( par.omega() ) > Constants::epsilon ) rc = 1. / fabs( par.omega() );
            double rw     = layer->rMid();
            double alpha  = acos( 1 - rw * rw / ( 2 * rc * rc ) );
            double fltLen = rw;
            if ( fabs( 1 - rw * rw / ( 2 * rc * rc ) ) < 1 ) fltLen = rc * alpha;
            z        = par.z0() + fltLen * par.tanDip();
            tof      = fltLen / Constants::c;
            double x = layer->getWire( corrWire )->xWireDC( z );
            double y = layer->getWire( corrWire )->yWireDC( z );
            delx     = -d0 * sinPhi0 - x;
            dely     = d0 * cosPhi0 - y;
            if ( m_tuplePick ) t_phiHit = thisHit->phi( z );
          }
          else
          {
            delx = -d0 * sinPhi0 - thisHit->x();
            dely = d0 * cosPhi0 - thisHit->y();
            if ( m_tuplePick ) t_phiHit = thisHit->phi();
          }
          predDrift = cosPhi0 * dely - sinPhi0 * delx + curv * ( delx * delx + dely * dely );
 
          // predDrift = predDrift * (1. - curv() * predDrift);
          ambig = ( predDrift >= 0 ) ? 1 : -1;
          if ( hasCurled ) ambig = -ambig;
          double entranceAngle = 0.;
          driftDist = thisHit->driftDist( tof + bunchTime, ambig, entranceAngle, 0., z );
          if ( alink == 0 )
          {
            // FIXME: is this ambig here VVVVV OK for incoming tracks?
            resid  = ambig * fabs( driftDist ) - predDrift;
            aresid = fabs( resid );
            // if (aresid > tkParam.maxInactiveResid ) ambig = 0;//yzhang delete 2012-10-09
          }
        }
        else
        {
          ambig = alink->ambig();
          if ( 6 == tkParam.lPrint ) cout << " pick ambig for hot" << endl;
          if ( m_tuplePick ) t_phiHit = par.phi0() + par.omega() * alink->fltLen();
        }
      }
 
      // yzhang 2012-08-20 , guess phi of this hit on z
      double   zGuess = par.z0() + layer->rMid() * par.tanDip();
      double   phiDCz = layer->getWire( corrWire )->phiDC( zGuess );
      BesAngle phiDCzMax( phiDCz + deltaPhiCellWidth );
      BesAngle phiDCzMin( phiDCz - deltaPhiCellWidth );
 
      if ( m_tuplePick && alink == 0 )
      {
        double sigma = 999.;
        if ( thisHit != 0 && alink == 0 )
        {
          double entranceAngle = 0.;
          sigma = thisHit->sigma( driftDist, ambig, entranceAngle, atan( par.tanDip() ), z );
        }
        m_pickPredDrift[t_iHit] = predDrift;
        m_pickDrift[t_iHit]     = driftDist;
        m_pickDriftTruth[t_iHit] =
            haveDigiDrift[thisHit->layernumber()][thisHit->wirenumber()];
        if ( curv * Bz <= 0. )
        {
          // positive track under minus Bz
          if ( ( phiDCzMin - phiExit > 0 ) || ( phiDCzMax - phiEnter < 0 ) )
            m_pickPhizOk[t_iHit] = 0;
          else m_pickPhizOk[t_iHit] = 1;
        }
        else
        {
          if ( ( phiDCzMin - phiEnter > 0 ) || ( phiDCzMax - phiExit < 0 ) )
            m_pickPhizOk[t_iHit] = 0;
          else m_pickPhizOk[t_iHit] = 1;
        }
        m_pickZ[t_iHit]     = z;
        m_pickWire[t_iHit]  = thisHit->wirenumber();
        m_pickResid[t_iHit] = aresid;
        if ( abs( sigma ) > 0.000001 ) m_pickSigma[t_iHit] = sigma;
        else m_pickSigma[t_iHit] = 999.;
        double t_phiLowCut  = -999.;
        double t_phiHighCut = -999.;
        if ( t_phiHit > -998. )
        {
          t_phiLowCut  = ( phiEnter - t_phiHit ) * rEnter;
          t_phiHighCut = ( phiExit - t_phiHit ) * rExit;
        }
        m_pickPhiLowCut[t_iHit]  = t_phiLowCut;
        m_pickPhiHighCut[t_iHit] = t_phiHighCut;
        m_pickDriftCut[t_iHit]   = goodDriftCut;
        m_pickMcTk[t_iHit]       = mcTkId;
        m_pickPt[t_iHit]         = tkFit->pt();
        m_pickCurv[t_iHit]       = curv;
        t_iHit++;
      }
 
      if ( curv * Bz <= 0. )
      {
        // positive track
        if ( ( phiDCzMin - phiExit > 0 ) || ( phiDCzMax - phiEnter < 0 ) )
        {
          if ( 6 == tkParam.lPrint )
          { std::cout << " CUT by phiDCz not in phi En Ex range, curv>=0" << std::endl; }
          continue;
        }
      }
      else
      {
        // negtive track
        if ( ( phiDCzMin - phiEnter > 0 ) || ( phiDCzMax - phiExit < 0 ) )
        {
          if ( 6 == tkParam.lPrint )
          { std::cout << " CUT by phiDCz not in phi En Ex range, curv<0" << std::endl; }
          continue;
        }
      }
 
      // Cases
      // (0) pred drift dist > goodDriftCut, drop Hit
      if ( ambig != 0 && fabs( predDrift ) > goodDriftCut )
      { // yzhang add 2012-08-17
        if ( 6 == tkParam.lPrint )
        { cout << " predDrift " << predDrift << ">goodDriftCut " << goodDriftCut << endl; }
        continue;
      }
 
      // (1) No good hit, but track near cell-edge => do nothing and continue
      if ( ambig == 0 && fabs( predDrift ) > goodDriftCut )
      { // yzhang 2009-11-03 add 3factor, 2011-05-30 from 3 to 2,2012-08-17 from 2 to 1
        if ( 6 == tkParam.lPrint )
        {
          cout << "   ambig==0 && |predDirft| " << fabs( predDrift ) << "> goodDriftCut "
               << goodDriftCut << endl;
          cout << " No good hit, but track near cell-edge " << endl;
        }
        continue;
      }
 
      // (2) Hit found:
      if ( ambig != 0 )
      {
        // yzhang changed 2009-10-19
        // if resid> maxActiveSimga * sigma do not add hits to track
        double entranceAngle = 0.;
        double sigma =
            thisHit->sigma( driftDist, ambig, entranceAngle, atan( par.tanDip() ), z );
        double factor = tkParam.pickHitFactor;
        // yzhang delete 2012-10-09
        // if(!is2d fabs(par.tanDip())<2) factor = (2-par.tanDip())/2 * factor;
        double residCut = tkParam.maxActiveSigma * factor * sigma; // yzhang 2012-08-21
        // if(6==tkParam.lPrint){
        // std::cout<< "aresid "<<aresid<<" residCut "<<residCut<<" sigma "<<sigma <<" tanDip
        // "<< par.tanDip() <<" factor "<<factor <<" tkParam.maxActiveSigma
        // "<<tkParam.maxActiveSigma<<std::endl;
        //}
 
        if ( alink == 0 && ( aresid <= residCut ) )
        {
          if ( 6 == tkParam.lPrint )
          {
            cout << " (2) New hit found " << endl; // yzhang debug
          }
          // yzhang 2012-08-17 delete
          int isActive = 1;
          // if (aresid > residCut) isActive = 0;
          // if(6==tkParam.lPrint) {
          //   if (aresid > residCut)
          //     std::cout << "notACT, resid  "<<aresid<<" >residCut " << residCut<< std::endl;
          // }
          MdcRecoHitOnTrack tempHot( *thisHit, ambig, bunchTime );
          tempHot.setActivity( isActive );
          // Don't add active hits if they're in use.
          if ( thisHit->usedHit() )
          {
            tempHot.setUsability( false );
            if ( 6 == tkParam.lPrint )
              std::cout << " thisHit used, setUsability false " << std::endl;
          }
          // very crude starting point for flight length !!!! improve
          double flt = layer->rMid();
          if ( hasCurled ) flt = Constants::twoPi * rCurl - layer->rMid();
          flt += 0.000001 * ( thisHit->x() + thisHit->y() );
          tempHot.setFltLen( flt );
          if ( 6 == tkParam.lPrint )
          {
            std::cout << " aresid  " << aresid << "<=residCut " << residCut << " nSig "
                      << aresid / sigma << " nSigCut " << ( tkParam.maxActiveSigma * factor )
                      << " factor " << factor << " maxActiveSigma " << tkParam.maxActiveSigma
                      << " tanDip " << par.tanDip() << std::endl;
            std::cout << " Append Hot  " << std::endl;
          }
          alink = (MdcHitOnTrack*)hitList->appendHot( &tempHot );
        }
        else
        {
          if ( 6 == tkParam.lPrint )
          {
            if ( alink != 0 ) { std::cout << "Exist hot found" << std::endl; }
            else if ( aresid > residCut )
            {
              thisHit->print( std::cout );
              std::cout << " Drop hit. aresid  " << aresid << ">residCut " << residCut
                        << " nSig " << aresid / sigma << " nSigCut "
                        << ( tkParam.maxActiveSigma * factor ) << " factor " << factor
                        << " maxActiveSigma " << tkParam.maxActiveSigma << " tanDip "
                        << par.tanDip() << std::endl;
            }
          }
        }
        if ( !localHistory.member( const_cast<MdcHitOnTrack*>( alink ) ) )
        {
          localHistory.append( const_cast<MdcHitOnTrack*>( alink ) );
          if ( hasCurled ) trk->setHasCurled();
          nFound++;
          if ( 6 == tkParam.lPrint )
            std::cout << " nFound=" << nFound << " nCand " << nCand << std::endl;
          if ( layer == firstInputLayer && firstInputHit < 0 ) { firstInputHit = nCand; }
        }
        else
        {
          if ( 6 == tkParam.lPrint )
            std::cout << "ErrMsg(warning) "
                      << "would have inserted identical HOT "
                         "twice in history buffer"
                      << std::endl;
        }
      }
 
      // (3) No hit found; if beyond known good region, should we continue?
      else if ( ambig == 0 && reachedLastInput )
      {
        if ( 6 == tkParam.lPrint ) std::cout << "ambig==0 " << std::endl;
        lContinue    = false;
        int nSuccess = 0;
        int last     = 8;
        if ( nCand < 8 ) last = nCand;
        for ( int i = 0; i < last; i++ )
        {
          int j = nCand - 1 - i;
          if ( localHistory[j] != 0 )
          {
            // std::cout<<  __FILE__ <<" localHistory["<<j<<"]!=0 nSuccess++ "<<  std::endl;
            nSuccess++;
          }
          if ( i == 2 && nSuccess >= 2 ) lContinue = true;
          if ( i == 4 && nSuccess >= 3 ) lContinue = true;
          if ( i == 7 && nSuccess >= 5 ) lContinue = true;
          if ( 6 == tkParam.lPrint )
            cout << i << "  (3) No hit found; if beyond known good region "
                 << endl; // yzhang debug
          if ( lContinue )
          {
            if ( 6 == tkParam.lPrint )
              std::cout << " pickHits break by lContinue. i=" << i << " nSuccess=" << nSuccess
                        << std::endl;
            break;
          }
        }
 
        if ( 6 == tkParam.lPrint )
          cout << "  (3) No hit found; if beyond known good region " << endl; // yzhang debug
        // if lContinue == false => there's been a gap, so quit
        if ( !lContinue )
        {
          // std::cout<<  __FILE__ <<"  "<< __LINE__  <<" break by !lContinue "<< std::endl;
          break;
        }
        localHistory.append( 0 );
      }
 
      nCand++;
      // Update last layer:
      if ( ambig != 0 && reachedLastInput )
      {
        if ( trk->hasCurled() == 0 )
        {
          if ( thisHit->layer()->rEnd() > trk->lastLayer()->rEnd() )
          { trk->setLastLayer( thisHit->layer() ); }
        }
        else
        {
          if ( thisHit->layer()->rEnd() < trk->lastLayer()->rEnd() )
          { trk->setLastLayer( thisHit->layer() ); }
        }
      }
 
    } // end loop over hit wires
    if ( t_iHit > 0 && m_tuplePick )
    {
      m_pickNCell         = tempDiff;
      m_pickNCellWithDigi = t_iHit;
      m_pickLayer         = layer->layNum();
      m_pickIs2D          = is2d;
      m_tuplePick->write();
    }
    if ( !lContinue && reachedLastInput )
    {
      // std::cout<<  __FILE__ <<" break by !lContinue && reachedLastInput "<< std::endl;
      break;
    }
  } // end loop over layers
 
  // Now look back at hits in early layer and see if there are any to be thrown away there.
  bool lContinue = true;
  for ( int ihit = firstInputHit; ihit >= 0; ihit-- )
  {
    if ( localHistory[ihit] != 0 )
    {
      if ( lContinue )
      {
        // Update firstLayer
        const MdcHitOnTrack* mdcHit = localHistory[ihit]->mdcHitOnTrack();
        if ( mdcHit != 0 )
        {
          if ( mdcHit->layer()->rEnd() < trk->firstLayer()->rEnd() )
          { trk->setFirstLayer( mdcHit->layer() ); }
        }
        continue; // no gap yet
      }
      else
      {
        // gap found; delete hits
        if ( 6 == tkParam.lPrint ) { std::cout << " gap found; delete hits. "; }
        if ( !localHistory[ihit]->isUsable() )
        {
          if ( 6 == tkParam.lPrint )
          {
            cout << "about to delete hit for unusable HOT:" << endl;
            localHistory[ihit]->print( std::cout );
          }
          hitList->removeHit( localHistory[ihit]->hit() );
        }
        if ( 6 == tkParam.lPrint )
        {
          cout << " current contents of localHistory: " << localHistory.length() << "hot"
               << endl;
          // for (int i=0; i<localHistory.length();++i) {
          // cout << " hit @ " << localHistory[i]->hit() << "; hot @ " << localHistory[i] <<
          // endl;
          // }
        }
        nFound--;
      }
    }
    else if ( localHistory[ihit] == 0 )
    {
      if ( 6 == tkParam.lPrint ) { cout << " localHistory= 0 " << endl; }
      int nSuccess = 0;
      lContinue    = false;
      int last     = 8;
      if ( nCand < 8 ) last = nCand;
      for ( int i = 0; i < last; i++ )
      {
        int j = ihit + 1 + i;
        if ( localHistory[j] != 0 ) nSuccess++;
        if ( i == 2 && nSuccess >= 2 ) lContinue = true;
        if ( i == 4 && nSuccess >= 3 ) lContinue = true;
        //  if (i == 7 && nSuccess >= 5) lContinue = true;
        if ( lContinue ) break;
      }
    }
  }
  if ( 6 == tkParam.lPrint )
  {
    std::cout << "nFound=" << nFound << " < " << tkParam.pickHitFract * nCand
              << " pickHitFract? " << tkParam.pickHitFract << "*" << nCand << std::endl;
  }
  if ( nFound < tkParam.pickHitFract * nCand )
    nFound = 0; // yzhang delete 2010-05-10 use pickHitFract=0.
  if ( 6 == tkParam.lPrint ) { cout << " localHistory= 0 " << endl; }
 
  if ( 6 == tkParam.lPrint || 3 == tkParam.lPrint )
  {
    cout << "After pickHits found " << nFound << " hit(s)" << endl;
    hitList->hotList().print( std::cout );
    std::cout << std::endl;
  }
 
  return nFound;
}
 
//************************************************************************
int MdcTrackList::finishHelix( MdcTrack& mdcTrk, const MdcHitMap* map,
                               const MdcDetector* gm ) {
  //***********************************************************************
  int success = 1;
 
  TrkRecoTrk&   trk = mdcTrk.track();
  TrkErrCode    fitResult;
  TrkHelixMaker helMaker;
  const TrkFit* tkFit = trk.fitResult();
  assert( tkFit != 0 );
  TrkHitList* hitList = trk.hits();
  assert( hitList != 0 );
  TrkExchangePar par = tkFit->helix( 0.0 );
 
  helMaker.setFlipAndDrop( trk, true, true );
  fitResult = hitList->fit();
 
  if ( !fitResult.success() && ( 3 == tkParam.lPrint ) )
  {
    cout << "Helix fit failure: " << endl;
    fitResult.print( cout );
  }
  helMaker.setFlipAndDrop( trk, false, false );
 
  if ( !fitResult.success() ) return 0;
 
  bool lcurler( fabs( tkFit->helix( 0 ).omega() ) > 3.4 ); // yzhang FIXME, ??
  pickHits( &mdcTrk, map, gm, lcurler );                   // yzhang add 2010-05-10
 
  if ( 3 == tkParam.lPrint )
    std::cout << __FILE__ << "   " << __LINE__ << " nHit after pickHit  " << hitList->nHit()
              << std::endl;
  if ( 3 == tkParam.lPrint ) hitList->hotList().printAll( std::cout );
  // if(hitList->nHit()<=0) return 0;
 
  helMaker.setFlipAndDrop( trk, true, true );
  fitResult = hitList->fit();
  if ( fitResult.failure() )
  {
    if ( 3 == tkParam.lPrint )
    {
      cout << "Second helix fit failed: " << endl;
      fitResult.print( std::cout );
    }
    return 0;
  }
  if ( 3 == tkParam.lPrint ) { cout << "Final fit: " << endl << trk << endl; }
  helMaker.setFlipAndDrop( trk, false, false );
 
  double chisqperDOF = 0.;
  int    nACT        = tkFit->nActive();
  int    nDOF        = nACT - 5;
  if ( nDOF > 5 ) { chisqperDOF = tkFit->chisq() / nDOF; }
  else { chisqperDOF = tkFit->chisq(); }
  if ( g_3dTkChi2 ) { g_3dTkChi2->fill( chisqperDOF ); } // yzhang hist cut
  if ( g_fitNAct ) { g_fitNAct->fill( nACT ); }          // yzhang hist cut
 
  int goodMatch = 1;
  if ( fitResult.failure() || chisqperDOF > tkParam.maxChisq || nACT < tkParam.minHits )
  {
    goodMatch = 0;
    if ( 3 == tkParam.lPrint )
    {
      std::cout << " goodMatch=0; chi2/dof " << chisqperDOF << " >?maxChisq"
                << tkParam.maxChisq << " nAct:" << nACT << "<?minHits" << tkParam.minHits
                << std::endl;
    }
  }
  // yzhang add
  ////yzhang FIXME
  // if (tkParam.lUseQualCuts) {
  // double tem2 = (float) trk.hits()->nHit() - nACT;
  // if (tem2 >= tkParam.maxNmissTrack) goodMatch = 0;
  // if (tem2 / float(trk.hits()->nHit()) > tkParam.maxNmissNorm)
  // goodMatch = 0;
  // }
  // zhangy add
 
  if ( goodMatch )
  {
    if ( 3 <= m_debug ) { std::cout << " ***** finishHelix success!" << std::endl; }
    trk.fitResult()->positionErr( 0.0 );
  }
  else
  { // Not a good match
    success = 0;
    if ( 3 <= m_debug ) { std::cout << " ***** finishHelix failure!" << std::endl; }
  } // end if goodmatch
 
  return success;
}
 
//************************************************************************
int MdcTrackList::finishCircle( MdcTrack& mdcTrk, const MdcHitMap* map,
                                const MdcDetector* gm ) {
  //************************************************************************
  TrkRecoTrk& trk = mdcTrk.track();
  if ( 9 == tkParam.lPrint )
  {
    std::cout << " finishCircle " << std::endl;
    trk.print( std::cout );
  }
 
  const TrkFit* tkFit = trk.fitResult();
  if ( 9 == tkParam.lPrint )
  { cout << "Before circle fit, nactive " << tkFit->nActive() << endl; }
  assert( tkFit != 0 );
  TrkHitList* hitList = trk.hits();
  assert( hitList != 0 );
  TrkCircleMaker circMaker;
  circMaker.setFlipAndDrop( trk, false, false ); // reset as a precaution
  // circMaker.setFactor(trk, 1.);//nSigma cut factor//yzhang FIXME 2010-09-13
 
  TrkErrCode fitResult = hitList->fit();
  if ( fitResult.failure() )
  {
    trk.status()->addHistory( TrkErrCode( TrkErrCode::fail, 15, "finishCircle" ),
                              "MdcTrkRecon" );
    if ( tkParam.lPrint > 1 )
    {
      cout << "First circle fit failed: " << endl;
      fitResult.print( std::cout );
    }
    return 0;
  }
  trk.status()->addHistory( TrkErrCode( TrkErrCode::succeed, 15, "finishCircle" ),
                            "MdcTrkRecon" );
 
  if ( 9 == tkParam.lPrint )
  { cout << "After circle fit, nactive " << tkFit->nActive() << endl; }
  double chisqperDOF;
  int    nDOF = tkFit->nActive() - 3;
  if ( nDOF > 3 ) { chisqperDOF = tkFit->chisq() / nDOF; }
  else { chisqperDOF = tkFit->chisq(); }
 
  if ( g_cirTkChi2 ) { g_cirTkChi2->fill( chisqperDOF ); } // yzhang hist cut
  int success = ( chisqperDOF <= tkParam.maxChisq && tkFit->nActive() >= 3 );
 
  if ( 9 == tkParam.lPrint )
  {
    std::cout << __FILE__ << " " << __LINE__ << " success " << success << " chisqperDOF "
              << chisqperDOF << "<? maxChisq " << tkParam.maxChisq << " nAct "
              << tkFit->nActive() << ">=3 " << std::endl;
    mdcTrk.track().hots()->printAll( std::cout );
  }
  bool lcurler( fabs( tkFit->helix( 0 ).omega() ) > 3.4 ); // yzhang FIXME, ??
  pickHits( &mdcTrk, map, gm, lcurler );
 
  if ( 9 == tkParam.lPrint )
  {
    std::cout << __FILE__ << "   " << __LINE__ << " nHit after pickHit  " << hitList->nHit()
              << std::endl;
  }
  // if(hitList->nHit()<=0) return 0;
 
  circMaker.setFlipAndDrop( trk, true, true );
  fitResult = hitList->fit();
  if ( fitResult.failure() )
  {
    if ( 9 == tkParam.lPrint )
    {
      cout << "Second circle fit failed: " << endl;
      fitResult.print( std::cout );
    }
    return 0;
  }
  if ( 9 == tkParam.lPrint ) { cout << "Final fit: " << endl << trk << endl; }
  circMaker.setFlipAndDrop( trk, false, false );
 
  nDOF = tkFit->nActive() - 3;
  if ( nDOF > 3 ) { chisqperDOF = tkFit->chisq() / nDOF; }
  else { chisqperDOF = tkFit->chisq(); }
  if ( g_cirTkChi2 ) { g_cirTkChi2->fill( chisqperDOF ); } // yzhang hist cut
  success = ( chisqperDOF <= tkParam.maxChisq ) && ( tkFit->nActive() >= 3 );
 
  if ( 9 == tkParam.lPrint )
  {
    cout << "chisqperDOF " << chisqperDOF << "=?" << tkParam.maxChisq << endl;
    cout << "nActive " << tkFit->nActive() << ">= 3" << endl;
  }
 
  if ( 9 == tkParam.lPrint ) { trk.printAll( cout ); }
 
  return success;
}
 
void MdcTrackList::dumpAxFill( const MdcTrack* trialTrack ) {
  std::cout << "ax fill: " << std::endl;
  if ( !trialTrack )
  {
    trialTrack->track().printAll( cout ); // yzhang debug
  }
}
 
void MdcTrackList::dumpSeed( const MdcSeg* seed, bool goodOnly ) {
  std::cout << std::endl << "Dump seed segment goodOnly=" << goodOnly << " ";
  seed->plotSegAll();
  std::cout << std::endl;
}
 
void MdcTrackList::dumpAxCombine( const MdcTrack* trialTrack ) {
  if ( NULL == trialTrack ) return;
  std::cout << std::endl << "-------------------------------------" << std::endl;
  std::cout << "Track and hitList after AxCombine " << std::endl;
  trialTrack->track().printAll( cout ); // yzhang debug
  TrkHotList::hot_iterator hotIter = trialTrack->track().hits()->hotList().begin();
  while ( hotIter != trialTrack->track().hits()->hotList().end() )
  {
    cout << "(" << ( (MdcHit*)( hotIter->hit() ) )->layernumber() << ","
         << ( (MdcHit*)( hotIter->hit() ) )->wirenumber() << ":" << hotIter->isActive()
         << ")  ";
    hotIter++;
  }
  std::cout << std::endl;
  std::cout << "-------------------------------------" << std::endl;
}
 
void MdcTrackList::dumpCircle( const MdcTrack* trialTrack ) {
  if ( NULL == trialTrack ) return;
  if ( !trialTrack->track().fitResult() ) return;
  /*
     double omega,r,pt;
     omega =trialTrack->track().fitResult()->helix(0.0).omega();
     if (omega == 0) pt = 0;
     else pt = (-1) * 1/(omega * 333.576 );
     std::cout<<"in MdcTrackList Circle Pt = "<< pt <<std::endl;//yzhang deubg
 
     if (omega == 0) r=0;
     else r = 1/ omega;
     std::cout<<"in MdcTrackList Circle R = "<< r <<std::endl;//yzhang deubg
     */
  std::cout << std::endl << "-------------------------------------" << std::endl;
  std::cout << "Track and hitList after finishCircle" << std::endl;
  trialTrack->track().printAll( cout );
  TrkHotList::hot_iterator hotIter = trialTrack->track().hits()->hotList().begin();
  while ( hotIter != trialTrack->track().hits()->hotList().end() )
  {
    cout << "(" << ( (MdcHit*)( hotIter->hit() ) )->layernumber() << ","
         << ( (MdcHit*)( hotIter->hit() ) )->wirenumber() << ":" << hotIter->isActive()
         << ")  ";
    hotIter++;
  }
  cout << endl;
  std::cout << "-------------------------------------" << std::endl;
}
 
void MdcTrackList::dumpD0( const TrkExchangePar& par ) {
  // yzhang hist
  // m_hd0->fill(par.d0());
  // m_d0 = par.d0();
  //   m_tuple1->write();//yzhang hist
  std::cout << std::endl << " Dump d0()  " << par.d0() << "\n"; // yzhang debug
 
  // zhangy hist
}
void MdcTrackList::dumpStFill() {
  std::cout << "Plot segs after st fillWithSegs " << std::endl;
  cout << endl;
}
 
void MdcTrackList::dumpStCombine( const MdcTrack* trialTrack ) {
  std::cout << std::endl << "-------------------------------------" << std::endl;
  std::cout << "Track and hitList after StCombine " << std::endl;
  if ( trialTrack ) trialTrack->track().printAll( cout );
  TrkHotList::hot_iterator hotIter = trialTrack->track().hits()->hotList().begin();
  int                      tmplay  = -1;
  while ( hotIter != trialTrack->track().hits()->hotList().end() )
  {
    int layer = ( (MdcHit*)( hotIter->hit() ) )->layernumber();
    if ( ( layer % 4 ) == 0 )
    {
      if ( tmplay != layer ) cout << endl;
    }
    cout << "(" << layer << "," << ( (MdcHit*)( hotIter->hit() ) )->wirenumber()
         << " act:" << hotIter->isActive() << " lr:" << hotIter->ambig() << ")  ";
    hotIter++;
    tmplay = layer;
  }
  cout << endl;
  std::cout << "-------------------------------------" << std::endl;
}
void MdcTrackList::dumpHelix( const MdcTrack* trialTrack ) {
  std::cout << std::endl << "-------------------------------------" << std::endl;
  std::cout << "Track and hitList after finishHelix " << std::endl;
  trialTrack->track().printAll( cout );
  TrkHotList::hot_iterator hotIter = trialTrack->track().hits()->hotList().begin();
  int                      tmplay  = -1;
  while ( hotIter != trialTrack->track().hits()->hotList().end() )
  {
    int layer = ( (MdcHit*)( hotIter->hit() ) )->layernumber();
    if ( ( layer % 4 ) == 0 )
    {
      if ( tmplay != layer ) cout << endl;
    }
    cout << "(" << layer << "," << ( (MdcHit*)( hotIter->hit() ) )->wirenumber() << ":"
         << hotIter->isActive() << ")  ";
    hotIter++;
    tmplay = layer;
  }
  cout << endl;
  std::cout << "-------------------------------------" << std::endl;
}
 
void MdcTrackList::dropMultiHotInLayer( const MdcTrack* tk ) {
  double tdr[43];
  double tdr_wire[43];
  for ( int i = 0; i < 43; i++ ) { tdr[i] = 9999.; }
 
  // make flag
  TrkHotList::hot_iterator hotIter = tk->track().hits()->hotList().begin();
  while ( hotIter != tk->track().hits()->hotList().end() )
  {
    MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*>( &( *hotIter->mdcHitOnTrack() ) );
 
    // driftTime(tof,z)
    double dt    = hot->mdcHit()->driftTime( 0., 0. );
    int    layer = hot->mdcHit()->layernumber();
    int    wire  = hot->mdcHit()->wirenumber();
    if ( dt < tdr[layer] )
    {
      tdr[layer]      = dt;
      tdr_wire[layer] = wire;
    }
    hotIter++;
  }
 
  std::cout << " tdr wire ";
  for ( int i = 0; i < 43; i++ ) { std::cout << i << " " << tdr[i] << " " << tdr_wire << " "; }
  std::cout << " " << std::endl;
  // inactive multi hit
  hotIter = tk->track().hits()->hotList().begin();
  while ( hotIter != tk->track().hits()->hotList().end() )
  {
    int    layer = hotIter->mdcHitOnTrack()->mdcHit()->layernumber();
    int    wire  = hotIter->mdcHitOnTrack()->mdcHit()->wirenumber();
    double dt    = hotIter->mdcHitOnTrack()->mdcHit()->driftTime( 0., 0. );
 
    if ( ( tdr[layer] < 9998. ) && ( tdr_wire[layer] != wire ) )
    {
      MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*>( &( *hotIter->mdcHitOnTrack() ) );
      hot->setActivity( false );
 
      std::cout << __FILE__ << " inactive " << layer << " " << wire << " dt " << dt
                << std::endl;
    }
    hotIter++;
  }
}