MdcxFindTracks.cxx

Go to the documentation of this file.

//-------------------------------------------------------------------------
// File and Version Information:
//  $Id: MdcxFindTracks.cxx,v 1.19 2011/12/08 06:52:29 zhangy Exp $
//
// Description:
//  Class Implementation for track finding
//
// Environment:
//  Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
//      S.  Wagner
//      Zhang Yao(zhangyao@ihep.ac.cn) Migration for BESIII
//
// Copyright Information:
//  Copyright (C) 1995  BEPCII
//
// History:
//  Migration for BESIII MDC
//
//------------------------------------------------------------------------
#include "MdcxReco/MdcxFindTracks.h"
#include "MdcGeom/Constants.h"
#include "MdcxReco/MdcxHit.h"
#include "MdcxReco/MdcxParameters.h"
#include "MdcxReco/Mdcxprobab.h"
#include <math.h>
 
#include "AIDA/IHistogram1D.h"
#include "AIDA/IHistogram2D.h"
#include "GaudiKernel/NTuple.h"
 
#define GET_NAME( n ) #n
#define prt( n ) << setw( 8 ) << GET_NAME( n ) << " " << setw( 8 ) << n << endl
 
extern int                 g_eventNo;
extern AIDA::IHistogram1D* g_omegag;
extern AIDA::IHistogram2D* g_dPhiAU;
extern AIDA::IHistogram2D* g_dPhiAV;
extern AIDA::IHistogram1D* g_trkllmk;
extern AIDA::IHistogram1D* g_trklcircle;
extern AIDA::IHistogram1D* g_trklgood;
extern AIDA::IHistogram1D* g_trklhelix;
extern AIDA::IHistogram1D* g_trkldrop1;
extern AIDA::IHistogram1D* g_trkldrop2;
extern AIDA::IHistogram1D* g_trklappend1;
extern AIDA::IHistogram1D* g_trklappend2;
extern AIDA::IHistogram1D* g_trklappend3;
extern AIDA::IHistogram1D* g_trklfirstProb;
extern AIDA::IHistogram1D* g_trkltemp;
extern AIDA::IHistogram1D* g_trklproca;
extern AIDA::IHistogram1D* g_trkld;
extern AIDA::IHistogram1D* g_trkle;
extern AIDA::IHistogram1D* g_trkldoca;
extern AIDA::IHistogram1D* g_trkllayer;
extern AIDA::IHistogram2D* g_trkldl;
extern AIDA::IHistogram2D* g_trklel;
extern AIDA::IHistogram2D* g_dropHitsSigma;
// extern AIDA::IHistogram1D*  g_trklprocaSl;
extern NTuple::Tuple*       m_xtupleAddSeg1;
extern NTuple::Item<long>   m_addSegSame;
extern NTuple::Item<double> m_addSegSeedSl;
extern NTuple::Item<double> m_addSegSeedPhi;
extern NTuple::Item<double> m_addSegSeedPhiLay;
extern NTuple::Item<double> m_addSegSeedD0;
extern NTuple::Item<double> m_addSegSeedLen;
extern NTuple::Item<double> m_addSegSeedPhi0;
extern NTuple::Item<double> m_addSegAddSl;
extern NTuple::Item<double> m_addSegAddPhi;
extern NTuple::Item<double> m_addSegAddPhiLay;
extern NTuple::Item<double> m_addSegAddD0;
extern NTuple::Item<double> m_addSegAddLen;
extern NTuple::Item<double> m_addSegAddPhi0;
 
extern NTuple::Tuple*       m_xtupleAddSeg2;
extern NTuple::Item<long>   m_addSegEvtNo;
extern NTuple::Item<double> m_addSegPoca;
extern NTuple::Item<long>   m_addSegSlayer;
extern NTuple::Item<double> m_addSegLen;
 
extern NTuple::Tuple*       m_xtupleSegComb;
extern NTuple::Item<double> m_segCombOmega;
extern NTuple::Item<double> m_segCombSameAU;
extern NTuple::Item<double> m_segCombSameUV;
extern NTuple::Item<double> m_segCombDLenAU;
extern NTuple::Item<double> m_segCombDLenUV;
extern NTuple::Item<double> m_segCombSlA;
extern NTuple::Item<double> m_segCombSlU;
extern NTuple::Item<double> m_segCombSlV;
extern NTuple::Item<double> m_segCombPhiA;
extern NTuple::Item<double> m_segCombPhiU;
extern NTuple::Item<double> m_segCombPhiV;
extern NTuple::Item<double> m_segCombPoca;
 
extern NTuple::Tuple*       m_xtupleDropHits;
extern NTuple::Item<long>   m_segDropHitsEvtNo;
extern NTuple::Item<long>   m_segDropHitsLayer;
extern NTuple::Item<long>   m_segDropHitsWire;
extern NTuple::Item<double> m_segDropHitsPull;
extern NTuple::Item<double> m_segDropHitsDoca;
extern NTuple::Item<double> m_segDropHitsSigma;
extern NTuple::Item<double> m_segDropHitsDrift;
extern NTuple::Item<double> m_segDropHitsMcTkId;
 
// IOS_USING_DECLARATION_MARKER - BesIII iostreams migration, do not touch this line!
using std::cout;
using std::endl;
using std::ostream;
 
MdcxFindTracks::MdcxFindTracks() {
  probmin = MdcxParameters::minTrkProb;
  curvmax = 1000000.0;
}
 
MdcxFindTracks::MdcxFindTracks( const HepAList<MdcxSeg>& segList, int debug ) {
  probmin = MdcxParameters::minTrkProb;
  curvmax = 1000000.0;
  m_debug = debug;
  process( segList );
}
 
MdcxFindTracks::~MdcxFindTracks() { KillList(); }
 
void MdcxFindTracks::process( const HepAList<MdcxSeg>& segList ) {
  static int   nfound  = 0;
  static float sumprob = 0.0;
 
  MdcxHel thel( 0., 0., 0., 0., 0. );
  mhel     = thel;
  int nSeg = segList.length();
  for ( int i = 0; i < nSeg; i++ ) { segList[i]->SetUsedOnHel( 0 ); }
 
  double xp, yp, xl1, yl1, rl1, xl2, yl2, rl2;
  double d0g, phi0g, omegag, z0g, tanlg;
  double d0g_sl, phi0g_sl, omegag_sl, z0g_sl, tanlg_sl;
  double zintAU, zintAV, zintAU_sl, zintAV_sl = 9999.;
  double rl1_sl, rl2_sl;
  double xc = 1.e9, yc = 1.e9;// add initialize 25-05-15
  double sinPhi0g_sl, cosPhi0g_sl, xp_sl, yp_sl;
  double phiAx = 1.e9, phiStU, phiStV, xl1_sl, yl1_sl, xl2_sl, yl2_sl, xl1_0, yl1_0, xl2_0, yl2_0;// add initialize 25-05-15
  double sx1, sy1, sx2, sy2;
  double x0g, y0g;
  double fltAx = 1.e9, ztest, fltL_sl, ztest_sl;// add initialize 25-05-15
  int    floop, trkCount = 0;
 
  // The main loop; first get the 3 SL combo we'll use to make a trial helix
  for ( int iSegCombo = 0; iSegCombo < MdcxParameters::nSegCombo; iSegCombo++ )
  {
    // FIXME optimize: group segList by superlayer ahead
    int axlay  = MdcxParameters::findTrkGroup[iSegCombo][0];
    int stUlay = MdcxParameters::findTrkGroup[iSegCombo][1];
    int stVlay = MdcxParameters::findTrkGroup[iSegCombo][2];
 
    if ( 4 == m_debug )
      std::cout << " Test combo: (" << axlay << "," << stUlay << "," << stVlay << ")"
                << std::endl;
 
    float max_dPhiAU = MdcxParameters::maxDp12[iSegCombo];
    float max_dPhiAV = MdcxParameters::maxDp13[iSegCombo];
 
    //---------------------------------------------------------------------
    //  Axial SL segment furnishes initial guess for d0, phi0, omega
    //---------------------------------------------------------------------
    for ( int iAx = 0; iAx < nSeg; iAx++ )
    {
      bool b_useGood_stU = true; // yzhang add 2009-10-21 17:20:45
      bool b_useGood_stV = true; // yzhang add 2009-10-21 17:20:45
      if ( ( segList[iAx]->GetUsedOnHel() != 0 ) || ( segList[iAx]->Fail() != 0 ) ||
           ( segList[iAx]->SuperLayer( 0 ) != axlay ) || ( segList[iAx]->Origin() != -1 ) )
        continue;
 
      if ( 4 == m_debug )
      {
        std::cout << "---------1.Ax seg------ No." << iAx << std::endl;
        segList[iAx]->printSegAll();
      }
 
      // Skip low pt track
      omegag = segList[iAx]->Omega();
      if ( fabs( omegag ) > MdcxParameters::maxTrkOmega )
      { // FIXME 0.2
        if ( 4 == m_debug ) std::cout << "SKIP by maxTrkOmega " << fabs( omegag ) << std::endl;
        continue;
      }
      if ( 4 == m_debug ) std::cout << " iAx omegag = " << omegag << std::endl;
      if ( g_omegag ) g_omegag->fill( omegag );
 
      const HepAList<MdcxHit>& hitsAx = segList[iAx]->XHitList();
 
      // Interpret segment as straight line also
      d0g_sl    = segList[iAx]->D0_sl_approx();
      phi0g_sl  = segList[iAx]->Phi0_sl_approx();
      omegag_sl = 0.0;
 
      // save c-tors; do it locally
      sinPhi0g_sl = -sin( phi0g_sl );
      cosPhi0g_sl = cos( phi0g_sl );
      xp_sl       = segList[iAx]->Xref() + sinPhi0g_sl * d0g_sl;
      yp_sl       = segList[iAx]->Yref() + cosPhi0g_sl * d0g_sl;
      d0g_sl      = yp_sl * cosPhi0g_sl + xp_sl * sinPhi0g_sl;
      fltL_sl     = segList[iAx]->Xref() * cosPhi0g_sl - segList[iAx]->Yref() * sinPhi0g_sl;
      if ( 4 == m_debug )
      {
        std::cout << "--Ax :" prt( d0g_sl ) prt( phi0g_sl ) prt( omegag_sl ) prt( sinPhi0g_sl )
                         prt( xp_sl ) prt( yp_sl ) prt( fltL_sl )
                  << std::endl;
      }
 
      if ( fabs( omegag ) > Constants::epsilon )
      {
        MdcxHel ohel = TakeToOrigin( *segList[iAx] );
        omegag       = ohel.Omega();
        phi0g        = ohel.Phi0();
        d0g          = ohel.D0();
        xc           = ohel.Xc();
        yc           = ohel.Yc();
        x0g          = ohel.X0();
        y0g          = ohel.Y0();
        phiAx        = atan2( y0g - yc, x0g - xc );
        fltAx        = dlen( -2, phi0g, -1, segList[iAx]->Phi0(), omegag );
        // fltAx = dlen(phi0g,segList[iAx]->Phi0(), omegag); //yzhang TEMP 2011-10-17
        if ( 4 == m_debug )
        {
          std::cout << "--Ax :"
                    << " ohel ";
          ohel.print();
        }
      }
 
      //---------------------------------------------------------------------
      // First stereo SL segement and axial segment furnish one z intersection
      //---------------------------------------------------------------------
      // for (int iStU = 0; iStU < nSeg; iStU++ ) //yzhang delete 2009-10-21 18:13:50
      for ( int iStU = -1; true; )
      { // yzhang add
        if ( !b_useGood_stU && ( iStU >= ( nSeg - 1 ) ) ) break;
        if ( b_useGood_stU && ( iStU >= ( nSeg - 1 ) ) )
        {
          iStU          = -1;
          b_useGood_stU = false;
        }
        iStU++;
        if ( ( segList[iAx]->GetUsedOnHel() != 0 ) || ( segList[iStU]->GetUsedOnHel() != 0 ) ||
             ( segList[iStU]->Fail() != 0 ) || ( segList[iStU]->SuperLayer( 0 ) != stUlay ) ||
             ( segList[iStU]->Origin() != -1 )
             // yzhang add  2009-10-21 17:21:55
             || ( b_useGood_stU && ( segList[iStU]->Nhits() < 4 ) ) )
        { continue; }
 
        if ( 4 == m_debug )
        {
          std::cout << " Test combo: AUV " << axlay << "," << stUlay << "," << stVlay
                    << std::endl;
          std::cout << "---------2.St U seg ------No." << iStU << std::endl; // yzhang debug
          segList[iStU]->printSeg();
          std::cout << " omega of seg iStU = " << segList[iStU]->Omega() << std::endl;
        }
 
        // Test dPhi between Ax slayer and StU slayer
        const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
        double                   dPhiAU  = fabs( hitsAx[0]->phiMid() - hitsStU[0]->phiMid() );
        if ( dPhiAU > Constants::pi ) dPhiAU = Constants::twoPi - dPhiAU;
        if ( 4 == m_debug )
        {
          if ( dPhiAU > max_dPhiAU )
          { cout << "SKIP by dPhiAU " << dPhiAU << " > " << max_dPhiAU << endl; }
          else { cout << "KEEP by dPhiAU " << dPhiAU << " <= " << max_dPhiAU << endl; }
          std::cout << " phi mid Ax " << hitsAx[0]->phiMid() << " StU " << hitsStU[0]->phiMid()
                    << std::endl;
          std::cout << std::endl;
        }
        if ( g_dPhiAU ) g_dPhiAU->fill( dPhiAU, segList[iStU]->SuperLayer() );
 
        if ( dPhiAU > max_dPhiAU ) { continue; }
 
        // Calculate z by Ax and StU segs
        xl1_0 = segList[iStU]->Xline_bbrrf();
        yl1_0 = segList[iStU]->Yline_bbrrf();
        sx1   = segList[iStU]->Xline_slope();
        sy1   = segList[iStU]->Yline_slope();
        if ( 4 == m_debug )
          std::cout prt( xl1_0 ) prt( yl1_0 ) prt( sx1 ) prt( sy1 ) prt( omegag ) << std::endl;
        if ( fabs( omegag ) > Constants::epsilon )
        {
          zintAU = findz_cyl( iAx, iStU, segList );
          xl1    = xl1_0 + sx1 * zintAU;
          yl1    = yl1_0 + sy1 * zintAU;
          rl1    = sqrt( xl1 * xl1 + yl1 * yl1 );
          phiStU = atan2( yl1 - yc, xl1 - xc );
          if ( 4 == m_debug )
            cout prt( zintAU ) prt( xl1 ) prt( yl1 ) prt( rl1 ) prt( phiStU ) << endl;
        }
        else { zintAU = -9999.; }
 
        // Calculate z_sl by Ax and StU segs
        zintAU_sl = findz_sl( iAx, iStU, segList );
        xl1_sl    = xl1_0 + sx1 * zintAU_sl;
        yl1_sl    = yl1_0 + sy1 * zintAU_sl;
        rl1_sl    = sqrt( xl1_sl * xl1_sl + yl1_sl * yl1_sl );
        if ( 4 == m_debug )
          cout prt( zintAU_sl ) prt( xl1_sl ) prt( yl1_sl ) prt( rl1_sl ) << endl;
 
        // Cut by z and z_sl
        if ( ( zintAU < -MdcxParameters::maxMdcZLen ) &&
             ( zintAU_sl < -MdcxParameters::maxMdcZLen ) )
        {
          if ( 4 == m_debug )
            std::cout << " SKIP by zintAU < max MdcZLen " << MdcxParameters::maxMdcZLen
                      << std::endl;
          continue;
        }
        if ( ( zintAU > MdcxParameters::maxMdcZLen ) &&
             ( zintAU_sl > MdcxParameters::maxMdcZLen ) )
        {
          if ( 4 == m_debug )
            std::cout << " SKIP by zintAU > max MdcZLen " << MdcxParameters::maxMdcZLen
                      << std::endl;
          continue;
        }
 
        //---------------------------------------------------------------------
        // Next stereo SL segement and axial segment furnish other z intersection
        //---------------------------------------------------------------------
        for ( int iStV = -1; true; )
        {
          if ( !b_useGood_stV && ( iStV >= ( nSeg - 1 ) ) )
          {
            if ( 4 == m_debug ) std::cout << iStV << " no segments in V slayer " << std::endl;
            break;
          }
          if ( b_useGood_stV && ( iStV >= ( nSeg - 1 ) ) )
          {
            iStV          = -1;
            b_useGood_stV = false;
          }
          iStV++;
          if ( ( segList[iStV]->GetUsedOnHel() != 0 ) ||
               ( segList[iStU]->GetUsedOnHel() != 0 ) ||
               ( segList[iAx]->GetUsedOnHel() != 0 ) || ( segList[iStV]->Fail() != 0 ) ||
               ( segList[iStV]->SuperLayer( 0 ) != stVlay ) ||
               ( segList[iStV]->Origin() != -1 )
               // yzhang add  2009-10-21 17:21:55
               || ( b_useGood_stV && ( segList[iStV]->Nhits() < 4 ) ) )
          { continue; }
 
          if ( 4 == m_debug )
          {
            std::cout << " Test combo: AUV " << axlay << "," << stUlay << "," << stVlay
                      << std::endl;
            std::cout << "---------3.St V seg ------No." << iStV << std::endl; // yzhang debug
            segList[iStV]->printSeg();
            std::cout << " omega of seg iStV = " << segList[iStV]->Omega() << std::endl;
          }
 
          // Calculate dPhi between Ax and StV segs
          const HepAList<MdcxHit>& hitsStV = segList[iStV]->XHitList();
          double                   dPhiAV = fabs( hitsAx[0]->phiMid() - hitsStV[0]->phiMid() );
          if ( dPhiAV > Constants::pi ) dPhiAV = Constants::twoPi - dPhiAV;
 
          if ( 4 == m_debug )
          {
            if ( dPhiAV > max_dPhiAV )
            { cout << "SKIP by dPhiAV " << dPhiAV << " > " << max_dPhiAV << endl; }
            else { cout << "KEEP by dPhiAV " << dPhiAV << " <= " << max_dPhiAV << endl; }
            std::cout << " phi mid Ax " << hitsAx[0]->phiMid() << " StV "
                      << hitsStV[0]->phiMid() << std::endl;
            std::cout << std::endl;
          }
          if ( g_dPhiAV ) g_dPhiAV->fill( dPhiAV, segList[iStV]->SuperLayer() );
          if ( dPhiAV > max_dPhiAV ) { continue; }
 
          // Calculate z by Ax and StV segs
          xl2_0 = segList[iStV]->Xline_bbrrf();
          yl2_0 = segList[iStV]->Yline_bbrrf();
          sx2   = segList[iStV]->Xline_slope();
          sy2   = segList[iStV]->Yline_slope();
          if ( fabs( omegag ) > Constants::epsilon )
          {
            zintAV = findz_cyl( iAx, iStV, segList );
            xl2    = xl2_0 + sx2 * zintAV;
            yl2    = yl2_0 + sy2 * zintAV;
            rl2    = sqrt( xl2 * xl2 + yl2 * yl2 );
            if ( 4 == m_debug )
            {
              segList[iAx]->printSeg();
              segList[iStV]->printSeg();
              cout << "zintAV " << zintAV << endl;
            }
            phiStV = atan2( yl2 - yc, xl2 - xc );
          }
          else { zintAV = -9999.; }
 
          // Calculate z by Ax and StV segs
          zintAV_sl = findz_sl( iAx, iStV, segList );
          xl2_sl    = xl2_0 + sx2 * zintAV_sl;
          yl2_sl    = yl2_0 + sy2 * zintAV_sl;
          rl2_sl    = sqrt( xl2_sl * xl2_sl + yl2_sl * yl2_sl );
 
          // Cut by z of AV segs
          if ( ( zintAV < -MdcxParameters::maxMdcZLen ) &&
               ( zintAV_sl < -MdcxParameters::maxMdcZLen ) )
          {
            if ( 4 == m_debug )
              std::cout << " SKIP by zintAV < max MdcZLen " << MdcxParameters::maxMdcZLen
                        << std::endl;
            continue;
          }
          if ( ( zintAV > MdcxParameters::maxMdcZLen ) &&
               ( zintAV_sl > MdcxParameters::maxMdcZLen ) )
          {
            if ( 4 == m_debug )
              std::cout << " SKIP by zintAV > max MdcZLen " << MdcxParameters::maxMdcZLen
                        << std::endl;
            continue;
          }
 
          // We now have enough info to form and fit a guess helix (inside iAx-iStU-iStV loops)
          if ( 4 == m_debug )
          {
            // cout << "KEEP layer[klm]: [" << axlay <<","<< stUlay ","<< stVlay << "] "
            // prt(d0g) prt(phi0g) prt(omegag) zrt(zintAU) prt(zintAV) << endl;
            segList[iAx]->printSeg();
            segList[iStU]->printSeg();
            segList[iStV]->printSeg();
          }
          MdcxFittedHel     fithel;
          double            helixFitSigma = MdcxParameters::helixFitSigma;
          double            first_prob    = 0.0;
          HepAList<MdcxHit> hitlist;
          hitlist.append( hitsAx );
          hitlist.append( hitsStU );
          hitlist.append( hitsStV );
          if ( g_trkllmk )
            for ( int i = 0; i < hitlist.length(); i++ )
            { g_trkllmk->fill( hitlist[i]->Layer() ); }
 
          //---------------------------------------------------------------------
          // Try fitting the trial helix from axial SL as circle first
          //---------------------------------------------------------------------
          floop = 1;
          while ( floop )
          {
            if ( 4 == m_debug )
              std::cout << "---------4.Ax circle fit------" << std::endl; // yzhang debug
            if ( fabs( omegag ) < Constants::epsilon )
            {
              if ( 4 == m_debug ) std::cout << "SKIP by omegag==0   " << std::endl;
              break;
            }
            if ( ( zintAU < -MdcxParameters::maxMdcZLen ) ||
                 ( zintAU > MdcxParameters::maxMdcZLen ) )
            {
              if ( 4 == m_debug )
                std::cout << "SKIP by zintAU out of range " << zintAU << " "
                          << MdcxParameters::maxMdcZLen << std::endl;
              break;
            }
            if ( ( zintAV < -MdcxParameters::maxMdcZLen ) ||
                 ( zintAV > MdcxParameters::maxMdcZLen ) )
            {
              if ( 4 == m_debug )
                std::cout << "SKIP by zintAU out of range " << zintAU << " "
                          << MdcxParameters::maxMdcZLen << std::endl;
              break;
            }
 
            // Calculate flight length of AUV segs
            if ( 4 == m_debug )
              cout << "dlen calc. slay U " << segList[iStU]->SuperLayer() << " slay V "
                   << segList[iStV]->SuperLayer() << endl;
            double fltLenUV =
                dlen( segList[iStU]->SuperLayer(), phiStU, segList[iStV]->SuperLayer(), phiStV,
                      omegag ); // yzhang change TEMP 2011-10-17
            // double fltLenUV = dlen(phiStU, phiStV, omegag);
            if ( 4 == m_debug )
            { cout prt( fltLenUV ) prt( phiStU ) prt( phiStV ) prt( omegag ) << std::endl; }
            if ( fltLenUV > MdcxParameters::maxDlen )
            { // 15. to 150.
              if ( 4 == m_debug )
                std::cout << "SKIP by fltLenUV > " << MdcxParameters::maxDlen << std::endl;
              break; // FIXME
            }
            tanlg = ( zintAV - zintAU ) / fltLenUV;
            if ( 4 == m_debug )
              cout << "dlen calc. slay A " << segList[iAx]->SuperLayer() << " slay U "
                   << segList[iStU]->SuperLayer() << endl;
            double fltLenAU =
                dlen( segList[iAx]->SuperLayer(), phiAx, segList[iStU]->SuperLayer(), phiStU,
                      omegag ); // yzhang change TEMP 2011-10-17
 
            if ( m_xtupleSegComb )
            {
              m_segCombOmega  = omegag;
              m_segCombSameAU = testFromSameTrack( segList[iAx], segList[iStU] );
              m_segCombSameUV = testFromSameTrack( segList[iStU], segList[iStV] );
              m_segCombSlA    = segList[iAx]->SuperLayer();
              m_segCombSlU    = segList[iStU]->SuperLayer();
              m_segCombSlV    = segList[iStV]->SuperLayer();
              m_segCombPhiA   = phiAx;
              m_segCombPhiU   = phiStU;
              m_segCombPhiV   = phiStV;
              m_segCombDLenUV = fltLenUV;
              m_segCombDLenAU = fltLenAU;
              m_xtupleSegComb->write();
            }
            // double fltLenAU = dlen(phiAx, phiStU, omegag);
            z0g = zintAU - fltLenAU * tanlg;
            if ( 4 == m_debug )
              cout prt( z0g ) prt( tanlg ) prt( fltLenAU ) prt( zintAU ) prt( zintAV ) << endl;
 
            // Cut by z0g
            if ( ( z0g < -MdcxParameters::maxMdcZLen ) ||
                 ( z0g > MdcxParameters::maxMdcZLen ) )
            {
              if ( 4 == m_debug )
                std::cout << "SKIP by z0g out of range " << z0g << ">"
                          << MdcxParameters::maxMdcZLen << std::endl;
              break;
            }
            ztest = z0g + fltAx * tanlg;
            if ( 4 == m_debug ) std::cout prt( ztest ) prt( fltAx ) << std::endl;
 
            //---Helix fitting
            // MdcxHel ghel(d0g,phi0g,omegag,z0g,tanlg); // ghel.print();
            MdcxHel ghel( segList[iAx]->D0(), segList[iAx]->Phi0(), segList[iAx]->Omega(),
                          ztest, tanlg, 0.0, 11111, 0, segList[iAx]->Xref(),
                          segList[iAx]->Yref() );
            // zoujh?: ghel.SetTurnFlag(1); // ghel.print();
            MdcxFittedHel firstfit( hitlist, ghel, helixFitSigma );
            first_prob = firstfit.Prob();
            if ( 4 == m_debug )
            {
              std::cout << " after first fit:  " << std::endl;
              firstfit.FitPrint();
            }
            if ( first_prob >= probmin )
            {
              MdcxHel       helixOrigin = TakeToOrigin( firstfit );
              MdcxFittedHel nextfit( hitlist, helixOrigin, helixFitSigma );
              first_prob = nextfit.Prob();
              fithel     = nextfit;
              if ( g_trklgood )
                for ( int i = 0; i < nextfit.Nhits(); i++ )
                { g_trklgood->fill( nextfit.Layer( i ) ); }
            }
            if ( g_trklfirstProb ) g_trklfirstProb->fill( first_prob );
            if ( 4 == m_debug ) cout << " prob after helix fitting = " << first_prob << endl;
            floop = 0;
          } // end fitting trial helix
 
          //---------------------------------------------------------------------
          // Try fitting straight line trial helix if fit from axial-SL-as-circle trial
          //  if helix not good enough
          //---------------------------------------------------------------------
          floop = 1;
          while ( floop )
          {
            if ( 4 == m_debug )
              std::cout << "---------4.2 straight line fit------" << std::endl; // yzhang debug
            if ( 4 == m_debug )
              cout << " zintAU_sl " << zintAU_sl << " zintAV_sl " << zintAV_sl << endl;
            if ( first_prob > probmin )
            {
              if ( 4 == m_debug ) cout << " helix fit good , exit straight line fit." << endl;
              break;
            }
 
            if ( ( zintAU_sl < -MdcxParameters::maxMdcZLen ) ||
                 ( zintAU_sl > MdcxParameters::maxMdcZLen ) )
              break;
            if ( ( zintAV_sl < -MdcxParameters::maxMdcZLen ) ||
                 ( zintAV_sl > MdcxParameters::maxMdcZLen ) )
              break;
 
            double dx = xl2_sl - xl1_sl;
            double dy = yl2_sl - yl1_sl;
            double dl = sqrt( dx * dx + dy * dy );
            tanlg_sl  = ( zintAV_sl - zintAU_sl ) / dl;
            dx        = xl1_sl + d0g_sl * sin( phi0g_sl );
            dy        = yl1_sl - d0g_sl * cos( phi0g_sl );
            dl        = sqrt( dx * dx + dy * dy );
            z0g_sl    = zintAU_sl - dl * tanlg_sl;
            if ( 4 == m_debug )
              cout prt( d0g_sl ) prt( phi0g_sl ) prt( z0g_sl ) prt( tanlg_sl ) << endl;
 
            if ( ( z0g_sl < -MdcxParameters::maxMdcZLen ) ||
                 ( z0g_sl > MdcxParameters::maxMdcZLen ) )
            {
              if ( 4 == m_debug )
                std::cout << "SKIP by z0g_sl:" << z0g_sl << " > " << MdcxParameters::maxMdcZLen
                          << std::endl;
              break;
            }
            // MdcxHel ghel(d0g_sl,phi0g_sl,omegag_sl,z0g_sl,tanlg_sl); // ghel.print();
            ztest_sl = z0g_sl + fltL_sl * tanlg_sl;
            if ( 4 == m_debug ) cout prt( ztest_sl ) << endl;
 
            // Helix fitting
            MdcxHel ghel( segList[iAx]->D0_sl_approx(), segList[iAx]->Phi0_sl_approx(),
                          omegag_sl, ztest_sl, tanlg_sl, 0.0, 11111, 0, segList[iAx]->Xref(),
                          segList[iAx]->Yref() );
            // zoujh?: ghel.SetTurnFlag(1); // ghel.print();
            MdcxFittedHel firstfit( hitlist, ghel, helixFitSigma );
            first_prob = firstfit.Prob();
            // if (first_prob >= probmin)fithel=firstfit;
            if ( 4 == m_debug )
            {
              ghel.print();
              std::cout << "first_prob " << first_prob << std::endl;
            }
 
            if ( first_prob >= probmin )
            {
              MdcxHel       helixOrigin = TakeToOrigin( firstfit );
              MdcxFittedHel nextfit( hitlist, helixOrigin, helixFitSigma );
              first_prob = nextfit.Prob();
              fithel     = nextfit;
              if ( 4 == m_debug )
              {
                cout << " three seg sl nextfit" << endl;
                nextfit.FitPrint();
              }
            }
            floop = 0;
          }
 
          //-----------------------------------------------------------------
          // Got a good trial helix fit; now try to add other segments onto it
          //-----------------------------------------------------------------
          floop = 1;
          while ( floop )
          {
            floop = 0;
            if ( 4 == m_debug )
              std::cout << "---------5. try add seg to helix------" << std::endl;
            if ( first_prob < probmin )
            {
              if ( 4 == m_debug )
                std::cout << "prob" << first_prob << " " << probmin
                          << " fit NOT good, exit add segs " << std::endl;
              break;
            }
            float bchisq = fithel.Chisq() * helixFitSigma * helixFitSigma;
            int   bndof  = fithel.Nhits() - 5;
            float bprob  = Mdcxprobab( bndof, bchisq );
            trkCount++;
            segList[iAx]->SetUsedOnHel( trkCount );
            segList[iStU]->SetUsedOnHel( trkCount );
            segList[iStV]->SetUsedOnHel( trkCount );
 
            if ( 4 == m_debug )
            {
              cout << "in add seg to trail helix, klm seg:" << endl; /// ZOUJH
              segList[iAx]->printSeg();                              /// ZOUJH
              segList[iStU]->printSeg();                             /// ZOUJH
              segList[iStV]->printSeg();                             /// ZOUJH
            }
 
            // Loop superlayers
            for ( int iSlay = 0; iSlay < Constants::nSuperLayer; iSlay++ )
            {
              int ilay = MdcxParameters::layerSet2AddSeg[iSegCombo][iSlay];
              for ( int iSeg = 0; iSeg < nSeg; iSeg++ )
              {
                if ( ( segList[iSeg]->SuperLayer( 0 ) != ilay ) ||
                     ( segList[iSeg]->GetUsedOnHel() != 0 ) ||
                     ( segList[iSeg]->Fail() != 0 ) || ( segList[iSeg]->Origin() != -1 ) )
                  continue;
                if ( 4 == m_debug )
                {
                  std::cout << endl << "== add seg " << iSeg << " ";
                  segList[iSeg]->printSeg();
                }
 
                // Cut by phi_diff = (Phi_Ax - Phi_Add)
                float  phiAxSeg     = segList[iAx]->XHitList()[0]->phiMid();
                float  phiAddSeg    = segList[iSeg]->XHitList()[0]->phiMid();
                double phiKNSegDiff = fabs( phiAxSeg - phiAddSeg );
                // yzhang add 2011-10-11
                // double phiKNSegDiff = (phiAxSeg - phiAddSeg);
                // if(phiKNSegDiff<-Constants::pi) phiKNSegDiff+=2*Constants::pi;
                // if(phiKNSegDiff>Constants::pi) phiKNSegDiff-=2*Constants::pi;
                int t_same = -999;
                if ( m_xtupleAddSeg1 )
                {
                  // std::cout<< __FILE__ << " testFromSameTrack  " << fromSameTrack<< "
                  // "<<std::endl;
                  t_same             = testFromSameTrack( segList[iAx], segList[iSeg] );
                  m_addSegSame       = t_same;
                  m_addSegSeedSl     = segList[iAx]->SuperLayer();
                  m_addSegSeedPhi    = segList[iAx]->XHitList()[0]->phiMid();
                  m_addSegSeedPhiLay = segList[iAx]->XHitList()[0]->Layer();
                  m_addSegSeedPhi0   = segList[iAx]->Phi0_sl_approx();
                  m_addSegSeedD0     = segList[iAx]->D0_sl_approx();
                  m_addSegAddSl      = segList[iSeg]->SuperLayer();
                  m_addSegAddPhi     = segList[iSeg]->XHitList()[0]->phiMid();
                  m_addSegAddPhiLay  = segList[iSeg]->XHitList()[0]->Layer();
                  m_addSegAddPhi0    = segList[iSeg]->Phi0_sl_approx();
                  m_addSegAddD0      = segList[iSeg]->D0_sl_approx();
                  m_xtupleAddSeg1->write();
                }
                // yzhang change 2011-10-11
                if ( phiKNSegDiff > 0.5 * Constants::pi && phiKNSegDiff < 1.5 * Constants::pi )
                // if (fabs(phiKNSegDiff) > 0.7*Constants::pi )
                { // yzhang TEMP
                  if ( 4 == m_debug )
                    std::cout << " SKIP by phi diff,same " << t_same << " Ax " << phiAxSeg
                              << " iSeg " << phiAddSeg << " diff " << phiKNSegDiff
                              << std::endl; // yzhang debug
                  continue;                 // zoujh FIXME
                }
                else
                {
                  if ( 4 == m_debug )
                    std::cout << " phi diff OK, Ax " << phiAxSeg << " added " << phiAddSeg
                              << " diff=" << phiKNSegDiff << std::endl; // yzhang debug
                }
 
                // Calculate doca
                xp = segList[iSeg]->Xref() - segList[iSeg]->SinPhi0() * segList[iSeg]->D0();
                yp = segList[iSeg]->Yref() + segList[iSeg]->CosPhi0() * segList[iSeg]->D0();
                const HepAList<MdcxHit>& hitsSegAdd = segList[iSeg]->XHitList();
                double proca = fithel.Doca( hitsSegAdd[0]->wx(), hitsSegAdd[0]->wy(),
                                            hitsSegAdd[0]->wz(), xp, yp );
                if ( g_trklproca ) g_trklproca->fill( proca );
 
                // Segment passes loose doca cut; does it belong on this track?
                if ( m_xtupleAddSeg2 )
                {
                  m_addSegEvtNo  = g_eventNo;
                  m_addSegPoca   = proca;
                  m_addSegSlayer = iSlay;
                  m_addSegLen    = fithel.Doca_Len();
                  m_xtupleAddSeg2->write();
                }
 
                if ( !( ( fabs( proca ) < MdcxParameters::maxProca ) &&
                        ( fithel.Doca_Len() < fithel.Lmax() ) ) )
                {
                  // if (g_trklprocaSl) g_trklprocaSl->fill(segList[iSeg]->SuperLayer(0));
                  if ( 4 == m_debug )
                  {
                    std::cout << " SKIP by fabs(proca):" << fabs( proca ) << "<"
                              << MdcxParameters::maxProca
                              << " or Doca_Len:" << fithel.Doca_Len() << "<" << fithel.Lmax()
                              << std::endl; // yzhang debug
                  }
                }
                else
                {
                  if ( 4 == m_debug )
                  {
                    std::cout << " proca & len OK, |proca| " << fabs( proca ) << "<"
                              << MdcxParameters::maxProca
                              << " && Doca_Len:" << fithel.Doca_Len() << "<" << fithel.Lmax()
                              << std::endl; // yzhang debug
                  }
                }
                if ( fithel.Doca_Len() < 0 ) continue; // yzhang add TEMP 2011-10-11
 
                if ( ( fabs( proca ) < MdcxParameters::maxProca ) &&
                     ( fithel.Doca_Len() < fithel.Lmax() ) )
                {
                  MdcxFittedHel newhel;
                  newhel.Grow( fithel, hitsSegAdd );
                  if ( newhel.Nhits() != hitlist.length() )
                  {
                    if ( 4 == m_debug )
                    {
                      cout << " rcs " << newhel.Rcs() << "<=?"
                           << MdcxParameters::maxRcsInAddSeg << endl;
                    }
                    // yzhang change 2009-10-21 FIXME
                    // if (newhel.Rcs() > MdcxParameters::maxRcsInAddSeg)
                    if ( newhel.Prob() < probmin )
                    {
                      if ( 4 == m_debug )
                      {
                        cout << " prob " << newhel.Prob() << "<" << probmin << ", ReFit"
                             << endl; /// ZOUJH
                      }
                      newhel.ReFit(); // FIXME
                    }
                    if ( 4 == m_debug )
                    {
                      cout << " refit prob " << newhel.Prob() << "<" << probmin << " Fail? "
                           << newhel.Fail() << endl;
                    }
                    // yzhang change 2009-10-21 FIXME
                    if ( ( newhel.Prob() >= probmin ) && ( newhel.Fail() == 0 ) )
                    {
                      // if ( (newhel.Rcs() <= MdcxParameters::maxRcsInAddSeg) &&
                      // (newhel.Fail() == 0) )
                      fithel  = newhel;
                      hitlist = newhel.XHitList();
                      segList[iSeg]->SetUsedOnHel( trkCount );
                      bchisq      = fithel.Chisq() * helixFitSigma * helixFitSigma;
                      bndof       = fithel.Nhits() - 5;
                      float tprob = Mdcxprobab( bndof, bchisq );
                      if ( tprob > bprob ) bprob = tprob;
                      if ( 4 == m_debug )
                      { cout << " segment ADDED, prob=" << newhel.Prob() << endl; }
                    } //*end trial track still good so keep new seg
                    else
                    {
                      if ( 4 == m_debug ) std::cout << " fit FAILED" << std::endl;
                    }
                  }
                  else
                  { //*end new hits so refit
                    segList[iSeg]->SetUsedOnHel( trkCount );
                    if ( 4 == m_debug ) cout << " segment ADDED, but no new hits" << endl;
                  } //*end dont need to refit if no new hits
                }   //*end fabs(proca) < 0.010
              }     //*end iSeg (loop over segs in a superlayer)
            }       //*end iSlay (loop over superlayers)
            if ( ( ( fithel.Prob() < 0.90 ) && ( fithel.Nhits() == 12 ) ) ||
                 ( fithel.Fail() != 0 ) )
            {
              for ( int i = 0; i < nSeg; i++ )
              {
                if ( segList[i]->GetUsedOnHel() == trkCount ) segList[i]->SetUsedOnHel( 0 );
              }
              trkCount--;
              break; //*jump to end floop first_prob >= probmin
            }
 
            if ( 4 == m_debug )
            {
              std::cout << "Track after add segs" << std::endl; // yzhang debug
              fithel.FitPrint();
            }
 
            // See if helix better moving backwards
            if ( 4 == m_debug )
              std::cout << "---------6. flip------" << std::endl; // yzhang debug
            fithel.flip();
 
            // See if this track shares lots of hits with another saved track
            int kki          = 0;
            int notduplicate = 1;
            while ( MdcxTrklist[kki] )
            {
              if ( 4 == m_debug )
                std::cout << "---------7. test saved track------" << std::endl; // yzhang debug
              const HepAList<MdcxHit>& junk  = MdcxTrklist[kki]->XHitList();
              int                      overl = 0;
              int                      kkj   = 0;
              while ( junk[kkj] )
              {
                int kkl = 0;
                while ( hitlist[kkl] )
                {
                  if ( hitlist[kkl] == junk[kkj] ) overl++;
                  kkl++;
                }
                kkj++;
              }
              if ( 4 == m_debug )
                cout << "overlap with track # " << kki << " is " << junk.length()
                     << " hitList " << hitlist.length() << " overlap " << overl << endl;
              if ( ( hitlist.length() / 4. ) <= overl ) notduplicate = 0;
              // yzhang 2009-10-21 16:30:53
              // if ( hitlist.length() <= (overl+2) ) notduplicate = 0;
              kki++;
            }
 
            if ( g_trklhelix )
              for ( int iHit = 0; iHit < fithel.Nhits(); iHit++ )
              { g_trklhelix->fill( fithel.Layer( iHit ) ); }
 
            // If not a duplicate, is it worth saving?
            if ( notduplicate )
            {
              if ( 4 == m_debug )
                std::cout << "---------8. test worth saving?------"
                          << std::endl; // yzhang debug
              MdcxFittedHel* finehel = new MdcxFittedHel( hitlist, fithel );
 
              if ( 4 == m_debug )
              {
                std::cout << __FILE__ << " finehel Prob>0.001 " << finehel->Prob() << " nhits "
                          << finehel->Nhits() << ">=15? "
                          << " bprob " << bprob << " >=?probmin " << probmin << " Fail==0? "
                          << finehel->Fail() << std::endl;
                finehel->FitPrint();
              }
              // cut nhits from 15->9 yzhang 2009-10-21
              if ( ( ( finehel->Prob() > 0.001 ) || ( finehel->Nhits() >= 15 ) ||
                     ( bprob > probmin )
                     //|| (bchisq/bndof > MdcxParameters::maxRcsInAddSeg)
                     ) &&
                   ( finehel->Fail() == 0 ) )
              {
                nfound++;
                sumprob += finehel->Prob();
                // if ( hitlist.length() == 16 )resout(finehel);
 
                // Begin hit dropping section
                // Set nodrop != 0 if want to prohibit dropping bad hits
                int nodrop = 0;
                if ( ( finehel->Prob() > 0.05 ) || nodrop )
                {
                  MdcxTrklist.append( finehel );
                  // yzhang hist
                  for ( int iHit = 0; iHit < finehel->Nhits(); iHit++ )
                  {
                    if ( g_trklappend1 )
                      g_trklappend1->fill( finehel->Layer( iHit ) ); // yzhang hist
                  }
                  // zhangy hist
                }
                else
                {
                  MdcxFittedHel* drophel = new MdcxFittedHel( drophits( finehel ) );
                  float ratdrop = float( drophel->Nhits() ) / float( finehel->Nhits() );
                  if ( 4 == m_debug )
                    cout << "APPEND track " << trkCount << ", drops "
                         << finehel->Nhits() - drophel->Nhits() << " helixnhit "
                         << finehel->Nhits() << " drophit " << drophel->Nhits()
                         << " hits, drop rate=" << ratdrop << ">?" << 0.5 << " prob "
                         << drophel->Prob() << " >?" << finehel->Prob() << " fail==0? "
                         << drophel->Fail() << endl;
                  if ( ( drophel->Prob() > finehel->Prob() ) && ( ratdrop > 0.50 ) &&
                       ( drophel->Fail() == 0 ) )
                  { // FIXME
                    if ( 4 == m_debug ) cout << "append drop helix " << endl;
                    MdcxTrklist.append( drophel );
                    if ( g_trklappend2 )
                      for ( int iHit = 0; iHit < drophel->Nhits(); iHit++ )
                      { g_trklappend2->fill( drophel->Layer( iHit ) ); }
                    delete finehel;
                  }
                  else
                  {
                    if ( 4 == m_debug ) cout << "append fine helix " << endl;
                    MdcxTrklist.append( finehel );
                    if ( g_trklappend3 )
                      for ( int iHit = 0; iHit < finehel->Nhits(); iHit++ )
                      { g_trklappend3->fill( finehel->Layer( iHit ) ); }
                    delete drophel;
                  }
                }
 
                // End hit dropping section
                int nwl = MdcxTrklist.length() - 1;
                if ( ( 4 == m_debug ) && ( nwl > -1 ) )
                {
                  if ( 4 == m_debug )
                  {
                    cout << endl << "found track " << trkCount << std::endl;
                    MdcxTrklist[nwl]->print();
                    MdcxTrklist[nwl]->FitPrint();
                  }
                }
              }
              else
              { //*endif good fit so add to MdcxTrklist
                for ( int i = 0; i < nSeg; i++ )
                {
                  if ( segList[i]->GetUsedOnHel() == trkCount ) segList[i]->SetUsedOnHel( 0 );
                }
                delete finehel;
                trkCount--;
              }
            } //*end if not duplicate
          }   //*end floop first_prob >= probmin
        }     // end iStV
      }       // end iStU
    }         // end iAx
    if ( 4 == m_debug ) cout << "end of this combo" << endl;
  } // end iSegCombo
  if ( 4 == m_debug )
    cout << " In MdcxFindTracks, found " << trkCount << " good tracks" << endl;
} // end of process
 
void MdcxFindTracks::print( ostream& o ) {
  o << "ngood:" << ngood << endl;
  o << "npure:" << npure << endl;
  o << "nbad:" << nbad << endl;
  o << "nfail:" << nfail << endl;
  o << "nhitsmc:" << nhitsmc << endl;
  o << "nhitsgd:" << nhitsgd << endl;
} // endof print
 
void MdcxFindTracks::plot() const {} // endof plot
 
void MdcxFindTracks::printpar( ostream& o ) {
  o << "|MdcxFindTracks| parameters" << endl;
  o << "probmin:" << probmin << endl;
  o << "curvmax:" << curvmax << endl;
} // endof printpar
 
void MdcxFindTracks::beginmore() {
  ngood   = 0;
  nbad    = 0;
  nfail   = 0;
  npure   = 0; // zero scalars
  nhitsmc = 0;
  nhitsgd = 0;      // zero scalars
  printpar( cout ); // print parameters to screen
} // endof beginmore
 
void MdcxFindTracks::endmore() { print( cout ); } // endof endmore
 
double MdcxFindTracks::findz_sl( int iAx, int iStU, const HepAList<MdcxSeg>& segList ) {
  double Ap = -sin( segList[iAx]->Phi0_sl_approx() );
  double Bp = cos( segList[iAx]->Phi0_sl_approx() );
  double xp = segList[iAx]->Xref() + Ap * segList[iAx]->D0_sl_approx();
  double yp = segList[iAx]->Yref() + Bp * segList[iAx]->D0_sl_approx();
  double xl = segList[iStU]->Xref() - segList[iStU]->SinPhi0() * segList[iStU]->D0();
  double yl = segList[iStU]->Yref() + segList[iStU]->CosPhi0() * segList[iStU]->D0();
 
  const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
  double                   Cl      = hitsStU[0]->wz();
  double                   ndotl   = Bp * hitsStU[0]->wy() + Ap * hitsStU[0]->wx();
  /*
     std::cout<<"ndotl   "<<ndotl<<" "<<Bp<<" "<<hitsStU[0]->wy()<<" "<<
     Ap<<" "<<hitsStU[0]->wx()<<std::endl;
     */
  double zint =
      ( ndotl != 0.0 ) ? ( ( Bp * ( yp - yl ) + Ap * ( xp - xl ) ) * Cl / ndotl ) : -1000.;
  // cout << " findz_sl segs " << iAx << " " << iStU << " zint " << zint << endl;
  return zint;
}
 
double MdcxFindTracks::findz_cyl( int iAx, int iStU, const HepAList<MdcxSeg>& segList ) {
  double zint = -1000.0;
  double r    = 1.0 / fabs( segList[iAx]->Omega() );
  double xl_0 = segList[iStU]->Xref() - segList[iStU]->SinPhi0() * segList[iStU]->D0();
  double yl_0 = segList[iStU]->Yref() + segList[iStU]->CosPhi0() * segList[iStU]->D0();
  double sx   = segList[iStU]->Xline_slope();
  double sy   = segList[iStU]->Yline_slope();
  double cx   = xl_0 - segList[iAx]->Xc();
  double cy   = yl_0 - segList[iAx]->Yc();
  double a    = sx * sx + sy * sy;
  double b    = 2.0 * ( sx * cx + sy * cy );
  // double c = (cx*cx + cy*cy - r*r);
  double bsq = b * b;
  double fac = 4.0 * a * ( cx * cx + cy * cy - r * r );
  double ta  = 2.0 * a;
  if ( fac < bsq )
  {
    double sfac = sqrt( bsq - fac );
    double z1   = -( b - sfac ) / ta;
    double z2   = -( b + sfac ) / ta;
    zint        = ( fabs( z1 ) < fabs( z2 ) ) ? z1 : z2;
  }
  // cout << " findz_cyl segs " << iAx << " " << iStU << " zint " << zint << endl;
  return zint;
}
 
double MdcxFindTracks::dlen( int slay1, double p1, int slay2, double p2, double om ) {
  // double dp = p2 - p1;
  double dp; // yzhang change 2011-10-17
  if ( ( slay1 == 2 || slay1 == 3 || slay1 == 4 || slay1 == 9 || slay1 == 10 ) )
  { // Ax
    dp = p2 - p1;
  }
  else if ( ( slay2 >= slay1 ) ) { dp = p2 - p1; }
  else { dp = p1 - p2; }
 
  if ( om < 0.0 )
  {
    while ( dp < -Constants::pi ) dp += Constants::pi;
    while ( dp > 0.0 ) dp -= 1 * Constants::pi;
  }
  else
  {
    while ( dp < 0.0 ) dp += 1 * Constants::pi;
    while ( dp > Constants::pi ) dp -= 1 * Constants::pi;
  }
  // double dl = dp * r;
  double dl = dp / om;
  // if(4 == m_debug) {
  //   double r = 1./om;
  //   cout prt(p1) prt(p2) prt(dp) prt(om) prt(r) prt(dl) << endl;
  // }
  return dl;
}
 
MdcxFittedHel MdcxFindTracks::drophits( MdcxFittedHel* finehel ) {
  MdcxHel                  mdcxHel   = (MdcxHel)*finehel;
  const HepAList<MdcxHit>& hitlist   = finehel->XHitList();
  HepAList<MdcxHit>        nwhitlist = hitlist;
  int                      nremove   = 0;
  int                      kkk       = 0;
  while ( nwhitlist[kkk] )
  {
    double pull  = nwhitlist[kkk]->pull( mdcxHel );
    int    layer = nwhitlist[kkk]->Layer();
 
    if ( m_xtupleDropHits )
    {
      float t_doca        = mdcxHel.Doca( *( nwhitlist[kkk] ) );
      float t_e           = nwhitlist[kkk]->e( t_doca );
      m_segDropHitsEvtNo  = g_eventNo;
      m_segDropHitsLayer  = layer;
      m_segDropHitsWire   = nwhitlist[kkk]->WireNo();
      m_segDropHitsPull   = pull;
      m_segDropHitsDoca   = t_doca;
      m_segDropHitsSigma  = t_e;
      m_segDropHitsDrift  = nwhitlist[kkk]->d( mdcxHel );
      m_segDropHitsMcTkId = nwhitlist[kkk]->getDigi()->getTrackIndex();
      m_xtupleDropHits->write();
    }
 
    if ( fabs( pull ) > MdcxParameters::dropHitsSigma[layer] )
    {
      nwhitlist.remove( kkk );
      nremove++;
    }
    else { kkk++; }
  }
  if ( m_debug == 4 )
    cout << " MdcxFindTracks::drophits drop " << nremove << " hits "
         << " nhits of drop helix = " << nwhitlist.length() << endl;
 
  MdcxFittedHel newhel( nwhitlist, mdcxHel );
  return newhel;
} // endof drophits
 
void MdcxFindTracks::resout( MdcxFittedHel* finehel ) {
  MdcxHel            mdcxHel = (MdcxHel)*finehel;
  HepAList<MdcxHit>& hitlist = (HepAList<MdcxHit>&)finehel->XHitList();
  int                kkk     = 0;
  while ( hitlist[kkk] )
  {
    hitlist[kkk]->SetConstErr( 0 );
    kkk++;
  }
  MdcxFittedHel newhel( hitlist, mdcxHel );
  if ( 4 == m_debug ) newhel.FitPrint( mhel, cout );
  kkk = 0;
  while ( hitlist[kkk] )
  {
    float doca = newhel.Doca( *hitlist[kkk] );
    float zh   = newhel.Doca_Zh();
    float tof  = newhel.Doca_Tof();
    float dd   = hitlist[kkk]->d( newhel );
    float tt   = hitlist[kkk]->tcor( zh, tof );
    if ( 4 == m_debug )
    {
      cout << "MdcxFindTracks::resout (" << hitlist[kkk]->Layer() << ","
           << hitlist[kkk]->WireNo() << ") dt " << tt << " resid " << fabs( dd ) - fabs( doca )
           << " doca " << fabs( doca ) << endl;
    }
    kkk++;
  }
  kkk = 0;
  while ( hitlist[kkk] )
  {
    hitlist[kkk]->SetConstErr( 1 );
    kkk++;
  }
} // endof drophits
 
MdcxHel MdcxFindTracks::TakeToOrigin( MdcxHel& ihel ) {
  double lng;
  double omegag = ihel.Omega();
  double phi0g  = ihel.Phi0();
  double phi0s  = phi0g;
  double xp     = ihel.Xref() - ihel.SinPhi0() * ihel.D0();
  double yp     = ihel.Yref() + ihel.CosPhi0() * ihel.D0();
  double d0g    = yp * cos( phi0g ) - xp * sin( phi0g );
  if ( omegag != 0.0 )
  {
    phi0g         = phi0g - Constants::pi;
    double xc     = sin( phi0g ) / omegag;
    double yc     = -cos( phi0g ) / omegag;
    double t1     = -xc - xp;
    double t2     = yc + yp;
    double phinew = atan2( t1, t2 );
    if ( omegag > 0.0 ) phinew += Constants::pi;
    if ( phinew > Constants::pi ) phinew -= Constants::twoPi;
    double xh  = xc - sin( phinew ) / omegag;
    double yh  = yc + cos( phinew ) / omegag;
    double x0  = xh + xp;
    double y0  = yh + yp;
    d0g        = sqrt( x0 * x0 + y0 * y0 );
    phi0g      = phinew + Constants::pi; ////???
    double sd0 = x0 * sin( phi0g ) - y0 * cos( phi0g );
    if ( sd0 > 0.0 ) d0g = -d0g;
    lng = dlen( -2, phi0g, -1, phi0s, omegag ); // yzhang TEMP 2011-10-17
    // lng = dlen(phi0g, phi0s, omegag);
  }
  else { lng = ihel.Xref() * ihel.CosPhi0() + ihel.Yref() * ihel.SinPhi0(); }
  double tanlg = ihel.Tanl();
  double z0g   = ihel.Z0() - lng * tanlg;
  // cout << " z0g, tanlg, lng " << z0g << " " << tanlg << " " << lng << endl;
  double  t0g   = ihel.T0();
  int     codeg = ihel.Code();
  MdcxHel ohel( d0g, phi0g, omegag, z0g, tanlg, t0g, codeg );
  return ohel;
}
 
bool MdcxFindTracks::testFromSameTrack( MdcxSeg* seg1, MdcxSeg* seg2 ) {
  /*
     std::cout<< " tkid seg1 "<<std::endl;
     for (int i =0; i<seg1->Nhits(); i++){
     int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
     std::cout<< tkid <<std::endl;
     }
     std::cout<< " tkid seg2 "<<std::endl;
     for (int i =0; i<seg2->Nhits(); i++){
     int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
     std::cout<< tkid <<std::endl;
     }
     */
 
  int trackId = -9999;
  for ( int i = 0; i < seg1->Nhits(); i++ )
  {
    int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
    if ( tkid >= 1000 ) tkid -= 1000;
    if ( ( trackId == -9999 ) )
    {
      if ( tkid >= 0 ) trackId = tkid;
    }
    else
    {
      if ( tkid != trackId ) return false;
    }
  }
  for ( int i = 0; i < seg2->Nhits(); i++ )
  {
    int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
    if ( tkid >= 1000 ) tkid -= 1000;
    if ( ( trackId == -9999 ) ) { return false; }
    else
    {
      if ( tkid != trackId ) return false;
    }
  }
  return true;
}