MdcTrackList Class Reference

#include <MdcTrackList.h>

Inheritance diagram for MdcTrackList:

List of all members.

Public Member Functions
	MdcTrackList (const MdcTrackParams &tkPar)
	~MdcTrackList ()
int	createFromSegs (MdcSegList *segs, const MdcHitMap *, const MdcDetector *, TrkContext &, double bunchTime)
int	finishCircle (MdcTrack &track, const MdcHitMap *, const MdcDetector *)
int	finishHelix (MdcTrack &track, const MdcHitMap *, const MdcDetector *)
int	pickHits (MdcTrack *, const MdcHitMap *, const MdcDetector *, bool pickAmb=true)
void	dumpSeed (const MdcSeg *seed, bool goodOnly)
void	dumpCircle (const MdcTrack *)
void	dumpAxFill (const MdcTrack *)
void	dumpAxCombine (const MdcTrack *)
void	dumpD0 (const TrkExchangePar &)
void	dumpStFill ()
void	dumpStCombine (const MdcTrack *)
void	dumpHelix (const MdcTrack *)
void	dropMultiHotInLayer (const MdcTrack *tk)
int	nTrack () const
void	setPlot (int plotFlag)
void	newParams (const MdcTrackParams &tkPar)
void	plot () const
void	store (RecMdcTrackCol *, RecMdcHitCol *)
int	arbitrateHits ()
void	remove (MdcTrack *atrack)
void	setD0Cut (double d0Cut)
void	setZ0Cut (double z0Cut)
void	setPtCut (double ptCut)
Static Public Attributes
static double	m_d0Cut = -999.
static double	m_z0Cut = -999.
static double	m_ptCut = -999.
Protected Attributes
MdcTrackParams	tkParam
Private Member Functions
MdcTrackList &	operator= (const MdcTrackList &)
	MdcTrackList (const MdcTrackList &)
Private Attributes
int	m_debug

---

Detailed Description

Definition at line 31 of file MdcTrackList.h.

---

Constructor & Destructor Documentation

MdcTrackList::MdcTrackList

(

const MdcTrackParams &

tkPar

)

Definition at line 91 of file MdcTrackList.cxx.

                                                       : 
  MdcTrackListBase(tkPar) { 
    // *************************************************************************
    return;
  }

MdcTrackList::~MdcTrackList

(

)

Definition at line 98 of file MdcTrackList.cxx.

{}

MdcTrackList::MdcTrackList

(

const MdcTrackList &

)

[private]

---

Member Function Documentation

int MdcTrackListBase::arbitrateHits

(

)

[inherited]

Definition at line 104 of file MdcTrackListBase.cxx.

References TrkHitList::begin(), TrkHitList::end(), TrkRecoTrk::fitResult(), for, MdcMap< K, V >::get(), TrkRecoTrk::hits(), genRecEmupikp::i, TrkRecoTrk::id(), MdcTrackParams::lPrint, MdcTrackParams::lRemoveInActive, MdcTrackListBase::nTrack(), MdcMap< K, V >::put(), q, TrkHitList::removeHit(), MdcTrackListBase::tkParam, and MdcTrack::track().

                                {
  //*************************************************************************
  // Look at all hits used in two or more tracks.  Assign hits to the track
  // that gives the lower residual.  If, however, many hits are shared by
  // a pair of tracks, assign them all to one or the other.
  // Refit any tracks that have had hits dropped.
  // The implementation is very clumsy, since the arrays were originally
  // indexed by id # => there is an unneeded layer of indexing.

  // return # of tracks deleted

  if (8 == tkParam.lPrint){
    std::cout << "=======Print before arbitrateHits=======" << std::endl;
  }

  int nDeleted = 0;
  std::vector<MdcTrack*> trksToKill;
  trksToKill.reserve(4);

  MdcMap<long,long> idMap;

  //usedInTrackNum records how many shared hits track has with each other track
  int* usedInTrackNum = new int [nTrack()];
  // to navigate from track id # to track pointer:
  MdcTrack** trkXRef = new MdcTrack* [nTrack()];
  //refitTrack flags track id #s of tracks to be refit
  int *refitTrack = new int [nTrack()];
  for (int i = 0; i < nTrack(); i++) {
    refitTrack[i] = 0;
  }

  // Fill xref table
  int itrack;
  for (itrack = 0; itrack < nTrack(); itrack++) {
    MdcTrack *atrack = (*this)[itrack];
    if (atrack == 0) continue;  // I don't think it can be, but . . .
    idMap.put(atrack->track().id(), itrack);
    trkXRef[itrack] = atrack;
  }
  // Loop through the tracks
  for (itrack = 0; itrack < nTrack(); itrack++) {

    if (8 == tkParam.lPrint) std::cout<<"arbitrate track No."<<itrack<<  std::endl;
    MdcTrack *atrack = (*this)[itrack];
    if (atrack == 0) continue;
    TrkRecoTrk& aRecoTrk = atrack->track();
    int lRefit = 0;
    int trackOld = -1;
    const TrkFit* tkFit = aRecoTrk.fitResult();
    assert (tkFit != 0);
    TrkHitList* hitList = aRecoTrk.hits();
    assert (hitList != 0);
restart:
    for (int ii = 0; ii < nTrack(); ii++) usedInTrackNum[ii] = 0;

    // Loop through hits on track, counting # used in other tracks
    int nPrev = 0;
    int nHitDeleted = 0;
    int maxGapLength = 0;//yzhang 2011-07-29 # of max continuous no hits layer for a track, Gap defined as missing layer >=2 
    int nGapGE2= 0;//yzhang 2011-07-29 # of no hits gap for a track 
    int nGapGE3= 0;//yzhang 2011-07-29 # of no hits gap for a track 
    int nHitInLayer[43];//yzhang 2010-09-20 for bad tracking testing
    int nDeleteInLayer[43];//yzhang 2010-09-20 
    for(int i=0;i<43;i++){
      nHitInLayer[i]=0;
      nDeleteInLayer[i]=0;
    }
    if(8 == tkParam.lPrint) std::cout<< "--arbitrate--"<<std::endl;
    for (TrkHitList::hot_iterator ihit(hitList->begin()); ihit != hitList->end(); ++ihit){
      int nUsed = ihit->hit()->nUsedHits();
      if (8 == tkParam.lPrint){
        std::cout<<"nUsed="<<nUsed<<":";
        ihit->hit()->printAll(std::cout);
      }
      if (8 == tkParam.lPrint) {
        double deltaChi = -999;
        ihit->getFitStuff(deltaChi);
        std::cout<< "deltaChi="<<deltaChi<<std::endl;
      }
      int layer = ihit->layerNumber();
      nHitInLayer[layer]++;

      if (!ihit->isActive()) {
        //-----------------------------------
        //yzhang delete not ACT hit 2010-05-14 
        //-----------------------------------
        if(tkParam.lRemoveInActive ) {//2010-05-16 
          nDeleteInLayer[layer]++;
          if (8 == tkParam.lPrint) {
            std::cout<< "=remove above inactive "<<std::endl;
          }
          TrkFundHit* hit = const_cast<TrkFundHit*> (ihit->hit());
          hitList->removeHit(hit);
          if(ihit == hitList->end()) break;
          --ihit;//be careful of the iterator, yzhang
        }
        continue;   // active hits only yzhang 2009-11-03 delete
      }
      if (nUsed > 1) {
        bool wasUsed = false;
        std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator> q =
          ihit->hit()->getUsedHits();
        for (TrkFundHit::hot_iterator i = q.first; i != q.second; ++i) {
          if ( !i->isActive() ) continue; //yzhang 2009-11-03 delete
          TrkRecoTrk * recoTrk=i->parentTrack();
          int id = recoTrk->id();
          if (id == aRecoTrk.id()) continue; //skip same track
          long index = 0;
          idMap.get(id, index);
          assert(index >= 0);
          usedInTrackNum[index]++;
          if (8 == tkParam.lPrint){
            std::cout<<" track "<<itrack<<"&" <<index
              << " shared hits "<<usedInTrackNum[index]<<":";
            ihit->printAll(std::cout);
          }
          wasUsed = true;
        }
        if (wasUsed) nPrev++;
      }// end nUsed > 1
    } // end loop over hits

    int testGap = 0;
    //std::cout<< __FILE__ << "   " << itrack<< "   "<<std::endl;
    for (int i=0;i<43;i++){
      //std::cout<< __FILE__ << "   " << i<< " nHitInLayer  "<<nHitInLayer[i]<<" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl;
      if (8 == tkParam.lPrint) {
        std::cout<<i<<" nHitInLayer  "<<nHitInLayer[i]
          <<" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl;
      }
      //1.only hit in layer deleted; 2.no hits in layer; 3.got hits in layer;
      if(nHitInLayer[i]>0 && (nHitInLayer[i]-nDeleteInLayer[i])==0) {
        //only hit in layer i has been deleted
        nHitDeleted++;
        if (8 == tkParam.lPrint) { 
          cout << "rec hits have been deleted in this layer"<<std::endl; 
        }
        testGap++;
        //std::cout<< __FILE__ << "   " << __LINE__ << " testGap3  "<<testGap<<std::endl;
      }else if(nHitInLayer[i]==0){
        //no hits in this layer i
        testGap++;
        //std::cout<< __FILE__ << "   " << __LINE__ << " testGap3  "<<testGap<<std::endl;
      }else{
        //std::cout<< __FILE__ << "   " << __LINE__ << " testGap3  "<<testGap<<std::endl;
        //got hit in layer i
        if(testGap>=2){
          nGapGE2++;
          if(testGap>=3){ nGapGE3++; }
          if(testGap>maxGapLength) maxGapLength=testGap;
          //std::cout<< __FILE__ << "   " << __LINE__ << " maxGapLength  "<<maxGapLength<<std::endl;
        }
        testGap=0;
      }//end for layer 43
    }

    bool toBeDeleted = false;

    if(tkParam.lPrint>1) std::cout<< "arbitrateHits tkNo:"<<itrack<<" nGapGE2=  "<<nGapGE2 << " nGapGE3=  "<<nGapGE3 << " maxGapLength=  "<<maxGapLength<<std::endl;
    //yzhang add nHitDeleted cut 2010-09-13 
    // remove track if # not Active 
    if (nHitDeleted >= tkParam.nHitDeleted) {
      if (tkParam.lPrint>1) {
        cout << "arbitrateHits: nHitDeleted "<<nHitDeleted<<" >= "<<tkParam.nHitDeleted
          <<" Killing tkNo " << itrack << endl;
      }
      toBeDeleted = true;
    }

    //yzhang add nGap cut 2011-07-29 
    // remove track with gaps and big gap
    if (nGapGE2 >= tkParam.nGapGE2) {
      if (tkParam.lPrint>1) {
        cout << "arbitrateHits: nGapGE2 "<<nGapGE2<<" >= "<<tkParam.nGapGE2 <<" Killing tkNo " << itrack << endl;
      }
      toBeDeleted = true;
    } 
    if (nGapGE3 >= tkParam.nGapGE3) {
      if (tkParam.lPrint>1) {
        cout << "arbitrateHits: nGapGE3 "<<nGapGE3<<" >= "<<tkParam.nGapGE3 <<" Killing tkNo " << itrack << endl;
      }
      toBeDeleted = true;
    } 
    if (maxGapLength >= tkParam.maxGapLength) {
      if (tkParam.lPrint>1) {
        cout << "arbitrateHits: maxGapLength "<<maxGapLength<<" >= "<<tkParam.maxGapLength<<" Killing tkNo " << itrack << endl;
      }
      toBeDeleted = true;
    } 

    if(toBeDeleted){
      nDeleted++;
      delete &(atrack->track());    // Delete the RecoTrk inside atrack
      atrack->setTrack(0);
      trksToKill.push_back(atrack);
      continue;
    } 

    //*******
    // How many hits are shared with a single track?
    int nMost = 0;
    int trackMost = 0;
    for (int ii = 0; ii < nTrack(); ii++) {
      if (8 == tkParam.lPrint){
        std::cout<<"tk:"<<itrack<<"&"<<ii
          <<" shared "<<usedInTrackNum[ii]<<" hits "<<  std::endl;
      }
      if (usedInTrackNum[ii] > nMost) {
        nMost = usedInTrackNum[ii];
        trackMost = ii;  //index of track w/ most hits in common w/ current trk
      }
    }

    // A little precaution against infinite loops:
    if (trackMost == trackOld) {
      std::cout << "ErrMsg(error) MdcTrackListBase:"
        << "Something ghastly happened in MdcTrackListBase::arbitrateHits"
        << std::endl;
      return 0;
    }
    trackOld = trackMost;


    //******
    // Decide whether to handle hits individually or in group
    double groupDiff = 0.0;   // relative quality of grouped hits for the two
    // tracks; > 0. => current track worse
    int nFound = 0;    // # of grouped hits located so far
    TrkHitOnTrk **theseHits = 0;  // grouped hits as seen in current track
    TrkHitOnTrk **thoseHits = 0;  // grouped hits as seen in the other track
    int lGroupHits = 0;

    if (nMost >= tkParam.nOverlap) {
      if (8 == tkParam.lPrint){
        std::cout<<"track "<<trackMost<<" shared "<<nMost<<" hits > Cut nOverlap "
          <<tkParam.nOverlap<<", group hits!"<<std::endl;
      }
      lGroupHits = 1;
      theseHits = new TrkHitOnTrk*[nMost];
      thoseHits = new TrkHitOnTrk*[nMost];
    }

    //*********
    // Go back through hits on this track, looking up the overlap of each
    // if grouping hits, only deal with hits shared with trackMost on this pass
    // otherwise, deal with all shared hits as encountered
    if(8 == tkParam.lPrint) std::cout<<"Go back through hits, looking up overlap hits"<< std::endl;
    if (nMost > 0) {
      if (8 == tkParam.lPrint) std::cout<<" nHits= "<< hitList->nHit()<< std::endl;
      for (TrkHitList::hot_iterator ihit(hitList->begin()); ihit != hitList->end(); ++ihit) {
        int nUsed = ihit->hit()->nUsedHits();

        if (8 == tkParam.lPrint){
          std::cout<< "--hit go back, nUsed="<<nUsed<<":";
          ihit->hit()->printAll(std::cout);
        }

        // only shared hits 
        if (nUsed < 2) { continue; }

        // active hits only
        if (!ihit->isActive()) {
          if (8 == tkParam.lPrint){ std::cout<<"act=0 continue"<<std::endl; }
          continue;   
        }

        //*** look at all overlaps for this hit
        std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator> q = ihit->hit()->getUsedHits();
        while (q.first!=q.second) { // nUsed > 0
          int dropThisHit = 0;
          TrkHitOnTrk *otherHot = const_cast<TrkHitOnTrk*>((--q.second).get());
          TrkRecoTrk *otherTrack = otherHot->parentTrack();

          if (!otherHot->isActive()) continue;

          // Again, skip "overlap" of track with itself
          if ( &aRecoTrk == otherTrack) continue;
          int otherId = otherTrack->id();
          long otherIndex = -1;
          idMap.get(otherId, otherIndex); assert(otherIndex >= 0);

          // if grouping hits, only look at hits shared with trackMost
          if (lGroupHits && otherIndex != trackMost) continue;

          if (lGroupHits) {
            if (8 == tkParam.lPrint) {
              std::cout<<"group hits "<<  std::endl;
            }
            // Calculate contribution of group to each chisq/dof
            //        groupDiff += fabs(ihit->resid(0)) -
            //          fabs(otherHot->resid(0));
            // Hack to handle tracks with 5 active hits:
            int aDof = tkFit->nActive() - 5;
            assert (otherTrack->fitResult() != 0);
            int otherDof = otherTrack->fitResult()->nActive() - 5;
            if (aDof <= 0) {groupDiff = 999;}
            else if (otherDof <= 0) {groupDiff = -999;}
            else {
              groupDiff += ihit->resid(0) * ihit->resid(0) * ihit->weight() /
                aDof -
                otherHot->resid(0) * otherHot->resid(0) * otherHot->weight() /
                otherDof;
            }
            theseHits[nFound] = const_cast<TrkHitOnTrk*>(ihit.get());
            thoseHits[nFound] = otherHot;
            nFound++;
            dropThisHit = 1;
          } else {   // handle hits individually

            if (8 == tkParam.lPrint) {
              std::cout<<"handle hits individually"<<  std::endl;
            }
            nFound++;
            if (fabs(ihit->resid(0)) > fabs(otherHot->resid(0)) ) {
              // turn off (inactivate) hit on this track
              lRefit = 1;
              //              ihit->hit()->setUnusedHit(ihit.get());
              //Should I be setting inactive, or deleting the hit???????
              const_cast<TrkHitOnTrk*>(ihit.get())->setActivity(0);
              dropThisHit = 1;
              if (8 == tkParam.lPrint) {
                std::cout<<"dorp hit ";
                const_cast<TrkHitOnTrk*>(ihit.get())->print(std::cout);
              }
              break;      // found other hit, so quit loop
            } else {
              // inactivate hit on other track
              refitTrack[otherIndex] = 1;
              //              otherHot->hit()->setUnusedHit(otherHot);
              otherHot->setActivity(0);
              if (8 == tkParam.lPrint) {
                std::cout<<"inactive hit on other track";
                const_cast<TrkHitOnTrk*>(ihit.get())->print(std::cout);
              }
              break;      // found other hit, so quit loop
            }
          } // end grouped/individual treatment

          if (dropThisHit == 1) break; // don't look for other matches since
          // this hit is now turned off
        } // end loop over nUsed

        // Quit if we've found all of the shared hits on this track
        if (lGroupHits && nFound == nMost  ||  nFound == nPrev) {
          if (8 == tkParam.lPrint) {
            std::cout<<"we've found all of the shared hits on this track,Quit"<<std::endl;
          }
          break;
        }

      }  // end loop over hits

      // Decide which track grouped hits belong with and inactivate accordingly
      if (lGroupHits) {
        if (8 == tkParam.lPrint) {
          cout << "nGroup: " << nMost << "  groupDiff: " << groupDiff << endl;
          cout << "Track: " << aRecoTrk.id() << "  nHit: "
            << hitList->nHit() << "  nActive: "
            << tkFit->nActive() << " chisq/dof: " <<
            tkFit->chisq()/(tkFit->nActive() - 5) << endl;
          TrkRecoTrk& othTrack = trkXRef[trackMost]->track();
          cout << "Track: "<< othTrack.id() << "  nHit: " <<
            othTrack.hits()->nHit() << "  nActive: " <<
            othTrack.fitResult()->nActive() << " chisq/dof: " <<
            othTrack.fitResult()->chisq() /
            (othTrack.fitResult()->nActive() - 5) << endl;
        }

        if (groupDiff > 0.0) {
          // inactivate hits on this track
          lRefit = 1;
          for (int ii = 0; ii < nMost; ii++) {
            TrkHitOnTrk *alink = theseHits[ii];
            TrkFundHit* hit = const_cast<TrkFundHit*> (alink->hit());
            hitList->removeHit(hit);//yzhang 2011-02-12 
            //alink->setActivity(0);
          }
          if (8 == tkParam.lPrint) std::cout<<"inactive hits on this track, No."<<aRecoTrk.id()<<  std::endl;
        } else {
          // inactivate hits on other track
          refitTrack[trackMost] = 1;
          for (int ii = 0; ii < nMost; ii++) {
            TrkHitOnTrk *alink = thoseHits[ii];
            TrkFundHit* hit = const_cast<TrkFundHit*> (alink->hit());
            hitList->removeHit(hit);//yzhang 2011-02-12 
            //alink->setActivity(0);
          }
          if (8 == tkParam.lPrint) std::cout<<"inactive hits on other track "<<   std::endl;
        }
        delete [] theseHits;
        delete [] thoseHits;

      } // end if lGroupHits

    } // end if nMost > 0

    //*********
    // Refit this track, if any hits have been dropped
    TrkErrCode fitResult;
    long index = -1;
    idMap.get(aRecoTrk.id(), index); assert (index >= 0);

    if (lRefit || refitTrack[index] == 1) {
      if (8 == tkParam.lPrint) {
        std::cout<<"after group ,refit track"<<aRecoTrk.id()<<  std::endl;
      }
      fitResult = hitList->fit();
      aRecoTrk.status()->addHistory(
          TrkErrCode(fitResult.success()?TrkErrCode::succeed:TrkErrCode::fail,14,"Arbitrated"), "MdcTrkRecon");
      if (fitResult.failure() && (8 == tkParam.lPrint )) {
        fitResult.print(std::cerr);
      }


      double chisqperDOF;
      bool badFit = true;
      if (fitResult.success()) {
        badFit = false;
        int nDOF = tkFit->nActive() - 5;
        if (nDOF > 5){
          chisqperDOF = tkFit->chisq() / nDOF;
        }else{
          chisqperDOF = tkFit->chisq();
        }

        if (chisqperDOF > tkParam.maxChisq) badFit = true;
        if (tkFit->nActive() < tkParam.minHits) badFit = true;
        double tem2 = (float) hitList->nHit() - tkFit->nActive();
        if (tkParam.lUseQualCuts) {
          if (tem2 >= tkParam.maxNmissTrack) badFit = true;
          if (tem2 /float(hitList->nHit()) > tkParam.maxNmissNorm){
            badFit = true;
          }
        }
        if(8== tkParam.lPrint) std::cout<<"fit quality:"<<
          " chisqperDof "<<chisqperDOF<<"?>"<<tkParam.maxChisq<<
            " nActive "<<tkFit->nActive()<<"?<"<<tkParam.minHits<<
            " nHit "<<hitList->nHit()<<" nhit-act "<<tem2<<"?>= nMiss "<<tkParam.maxNmissTrack<<
            " hit-act/nhit "<<tem2/float(hitList->nHit())<<"?> MissNorm "<<tkParam.maxNmissNorm
            <<  std::endl;


      }
      if (8 == tkParam.lPrint)  {
        cout << "Refitting track " << aRecoTrk.id() << " success = "
          << fitResult.success() << "\n";
      }
      // If the track no longer passes cuts, delete it
      if (fitResult.failure() || badFit ) {
        nDeleted++;
        //  Don't change the track list while we're iterating through it!
        //      remove(atrack);
        //int id = aRecoTrk.id();
        if (8 == tkParam.lPrint) {
          cout << "fitResult.failure? "<<fitResult.failure()
            <<" badFit? "<<badFit <<" Killing tkNo " << itrack << endl;
        }
        delete &(atrack->track());    // Delete the RecoTrk inside atrack
        atrack->setTrack(0);
        trksToKill.push_back(atrack);
        continue;
      }
    } // end if lRefit

    if (lGroupHits) goto restart;

  } // end loop over tracks
  if (8 == tkParam.lPrint) std::cout<<"end of loop over tracks"<<  std::endl;

  // Remove dead track husks
  for (int itk = 0; itk < (int)trksToKill.size(); itk++) {
    remove(trksToKill[itk]);
    if (8 == tkParam.lPrint) std::cout<<"remode dead track No."<<itk<<  std::endl;
  }
  if (8 == tkParam.lPrint) std::cout<<"---end of arbitrateHits"<<  std::endl;

  delete [] usedInTrackNum;
  delete [] refitTrack;
  delete [] trkXRef;
  return nDeleted;
}

int MdcTrackList::createFromSegs	(	MdcSegList *	segs,
		const MdcHitMap *	,
		const MdcDetector *	,
		TrkContext &	,
		double	bunchTime
	)			[virtual]

Implements MdcTrackListBase.

Definition at line 103 of file MdcTrackList.cxx.

References EvtCyclic3::append(), MdcTrackParams::combineByFitHits, MdcSegGrouper::combineSegs(), MdcSegList::count(), TrkExchangePar::d0(), dumpAxCombine(), dumpCircle(), dumpHelix(), dumpSeed(), MdcSegGrouper::dumpSegList(), dumpStCombine(), Constants::epsilon, MdcSegGrouperSt::fillWithSegs(), MdcSegGrouperAx::fillWithSegs(), finishCircle(), finishHelix(), TrkRecoTrk::fitResult(), MdcSegList::getSeed(), TrkFit::helix(), MdcTrackParams::lPrint, MdcTrackListBase::m_d0Cut, m_debug, MdcTrackListBase::m_ptCut, MdcTrackListBase::m_z0Cut, MdcTrackParams::maxSegChisqO, TrkAbsFit::pt(), MdcSegList::resetSeed(), MdcTrackListBase::tkParam, MdcTrack::track(), and TrkExchangePar::z0().

                                           {
  // *************************************************************************
  //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);

  // *** 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;

}

void MdcTrackList::dropMultiHotInLayer

(

const MdcTrack *

tk

)

Reimplemented from MdcTrackListBase.

Definition at line 1151 of file MdcTrackList.cxx.

References TrkHotList::begin(), dt, TrkHotList::end(), TrkRecoTrk::hits(), TrkHitList::hotList(), genRecEmupikp::i, TrkHitOnTrk::setActivity(), and MdcTrack::track().

                                                        {
  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++;
  }
}

void MdcTrackList::dumpAxCombine

(

const MdcTrack *

)

Definition at line 1057 of file MdcTrackList.cxx.

References TrkHotList::begin(), TrkHotList::end(), TrkRecoTrk::hits(), TrkHitList::hotList(), TrkRecoTrk::printAll(), and MdcTrack::track().

Referenced by createFromSegs().

                                                            {
  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::dumpAxFill

(

const MdcTrack *

)

Definition at line 1044 of file MdcTrackList.cxx.

References TrkRecoTrk::printAll(), and MdcTrack::track().

                                                          {
  std::cout << "ax fill: "<<std::endl; 
  if(!trialTrack){
    trialTrack->track().printAll(cout);//yzhang debug
  }
}

void MdcTrackList::dumpCircle

(

const MdcTrack *

)

Definition at line 1073 of file MdcTrackList.cxx.

References TrkHotList::begin(), TrkHotList::end(), TrkRecoTrk::fitResult(), TrkRecoTrk::hits(), TrkHitList::hotList(), TrkRecoTrk::printAll(), and MdcTrack::track().

Referenced by createFromSegs().

                                                        {
  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 &

)

Definition at line 1101 of file MdcTrackList.cxx.

References TrkExchangePar::d0().

                                                    {
  //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::dumpHelix

(

const MdcTrack *

)

Definition at line 1133 of file MdcTrackList.cxx.

References TrkHotList::begin(), TrkHotList::end(), TrkRecoTrk::hits(), TrkHitList::hotList(), TrkRecoTrk::printAll(), and MdcTrack::track().

Referenced by createFromSegs().

                                                       { 
  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::dumpSeed	(	const MdcSeg *	seed,
		bool	goodOnly
	)

Definition at line 1051 of file MdcTrackList.cxx.

References MdcSeg::plotSegAll().

Referenced by createFromSegs().

                                                              {
  std::cout << std::endl<<"Dump seed segment goodOnly="<<goodOnly<<" ";
  seed->plotSegAll();
  std::cout<< std::endl;
}

void MdcTrackList::dumpStCombine

(

const MdcTrack *

)

Definition at line 1116 of file MdcTrackList.cxx.

References TrkHotList::begin(), TrkHotList::end(), TrkRecoTrk::hits(), TrkHitList::hotList(), TrkRecoTrk::printAll(), and MdcTrack::track().

Referenced by createFromSegs().

                                                            {
  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::dumpStFill

(

)

Definition at line 1110 of file MdcTrackList.cxx.

                              {
  std::cout << "Plot segs after st fillWithSegs " << std::endl;
  cout <<endl;

}

int MdcTrackList::finishCircle	(	MdcTrack &	track,
		const MdcHitMap *	,
		const MdcDetector *
	)

Definition at line 951 of file MdcTrackList.cxx.

References TrkFitStatus::addHistory(), TrkAbsFit::chisq(), TrkErrCode::fail, TrkErrCode::failure(), TrkHitList::fit(), TrkRecoTrk::fitResult(), g_cirTkChi2, TrkFit::helix(), TrkRecoTrk::hits(), TrkRecoTrk::hots(), MdcTrackParams::lPrint, MdcTrackParams::maxChisq, TrkFit::nActive(), TrkHitList::nHit(), TrkExchangePar::omega(), pickHits(), TrkErrCode::print(), TrkRecoTrk::print(), TrkRecoTrk::printAll(), TrkHotList::printAll(), TrkSimpleMaker< T >::setFlipAndDrop(), TrkRecoTrk::status(), TrkErrCode::succeed, MdcTrackListBase::tkParam, and MdcTrack::track().

Referenced by createFromSegs().

                     {
  //************************************************************************
  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;
}

int MdcTrackList::finishHelix	(	MdcTrack &	track,
		const MdcHitMap *	,
		const MdcDetector *
	)

Definition at line 863 of file MdcTrackList.cxx.

References TrkAbsFit::chisq(), TrkErrCode::failure(), TrkHitList::fit(), TrkRecoTrk::fitResult(), g_3dTkChi2, g_fitNAct, TrkFit::helix(), TrkRecoTrk::hits(), TrkHitList::hotList(), MdcTrackParams::lPrint, m_debug, MdcTrackParams::maxChisq, MdcTrackParams::minHits, TrkFit::nActive(), TrkHitList::nHit(), TrkExchangePar::omega(), pickHits(), TrkAbsFit::positionErr(), TrkErrCode::print(), TrkHotList::printAll(), TrkSimpleMaker< T >::setFlipAndDrop(), TrkErrCode::success(), MdcTrackListBase::tkParam, and MdcTrack::track().

Referenced by createFromSegs().

                     {
  //***********************************************************************
  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
  //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;
}

void MdcTrackListBase::newParams

(

const MdcTrackParams &

tkPar

)

[inherited]

Definition at line 588 of file MdcTrackListBase.cxx.

References MdcTrackListBase::tkParam.

                                                       {
  //**************************************************************************
  tkParam = tkPar;
}

int MdcTrackListBase::nTrack

(

)

const [inline, inherited]

Definition at line 39 of file MdcTrackListBase.h.

Referenced by MdcTrackListBase::arbitrateHits(), MdcTrackListBase::plot(), and MdcTrackListBase::store().

{return length();}

MdcTrackList& MdcTrackList::operator=

(

const MdcTrackList &

)

[private]

int MdcTrackList::pickHits	(	MdcTrack *	,
		const MdcHitMap *	,
		const MdcDetector *	,
		bool	pickAmb = true
	)

Definition at line 295 of file MdcTrackList.cxx.

References abs, alpha, MdcHitOnTrack::ambig(), TrkHitList::appendHot(), TrkHitList::begin(), TrkRecoTrk::bField(), BField::bFieldZ(), Constants::c, cos(), TrkExchangePar::d0(), BesAngle::deg(), MdcHit::digi(), MdcHit::driftDist(), TrkHitList::end(), Constants::epsilon, MdcDetector::firstLayer(), MdcTrack::firstLayer(), TrkRecoTrk::fitResult(), TrkHitOnTrk::fltLen(), g_pickHitWire, TrkFundHit::getHitOnTrack(), RawData::getTrackIndex(), MdcTrack::hasCurled(), haveDigiDrift, TrkFit::helix(), TrkRecoTrk::hits(), MdcHitMap::hitWire(), genRecEmupikp::i, if(), TrkFitStatus::is2d(), ganga-rec::j, MdcDetector::lastLayer(), MdcTrack::lastLayer(), MdcHit::layer(), MdcHit::layernumber(), MdcTrackParams::lPrint, m_pickCurv, m_pickDrift, m_pickDriftCut, m_pickDriftTruth, m_pickIs2D, m_pickLayer, m_pickMcTk, m_pickNCell, m_pickNCellWithDigi, m_pickPhiHighCut, m_pickPhiLowCut, m_pickPhizOk, m_pickPredDrift, m_pickPt, m_pickResid, m_pickSigma, m_pickWire, m_pickZ, m_tuplePick, MdcTrackParams::maxActiveSigma, Constants::maxCell, MdcTrackParams::maxInactiveResid, TrkHitOnTrk::mdcHitOnTrack(), mdcWrapWire(), MdcDetector::nextLayer(), TrkExchangePar::omega(), MdcHit::phi(), TrkExchangePar::phi0(), Constants::pi, MdcTrackParams::pickHitFactor, MdcTrackParams::pickHitFract, MdcTrackParams::pickHitMargin, MdcTrackParams::pickHitPhiFactor, MdcTrackParams::pickSkipExistLayer, MdcTrackParams::pickUitlLastLayer, BesAngle::posRad(), MdcDetector::prevLayer(), MdcHit::print(), MdcTrack::projectToR(), TrkAbsFit::pt(), Constants::radToDegrees, TrkHitList::removeHit(), MdcLayer::rEnd(), MdcLayer::rIn(), MdcLayer::rOut(), TrkHitOnTrk::setActivity(), MdcTrack::setFirstLayer(), TrkHitOnTrk::setFltLen(), MdcTrack::setHasCurled(), MdcTrack::setLastLayer(), TrkHitOnTrk::setUsability(), MdcHit::sigma(), sin(), TrkRecoTrk::status(), TrkExchangePar::tanDip(), MdcTrackListBase::tkParam, MdcTrack::track(), TrkRecoTrk::trackT0(), Constants::twoPi, TrkFundHit::usedHit(), MdcHit::wirenumber(), MdcHit::x(), x, MdcHit::y(), and TrkExchangePar::z0().

Referenced by finishCircle(), and finishHelix().

                                        {

  /***************************************************************************/

  // 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();   
  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;
}

void MdcTrackListBase::plot

(

)

const [inherited]

Definition at line 92 of file MdcTrackListBase.cxx.

References MdcTrackListBase::nTrack().

                             {
  //*************************************************************************
  std::cout<< "nTrack   "<<nTrack() << std::endl;//yzhang debug
  for (int itrack = 0; itrack < nTrack(); itrack++) {
    MdcTrack *atrack = (*this)[itrack];
    if (atrack == NULL) continue;
    atrack->track().printAll(cout);
  }
}

void MdcTrackListBase::remove

(

MdcTrack *

atrack

)

[inherited]

Reimplemented in MdcTrackListCsmc.

Definition at line 595 of file MdcTrackListBase.cxx.

                                           {
  //--------------------------------------------------------------------
  if (atrack != 0) {
    HepAList<MdcTrack>::remove( atrack );
    delete atrack;
  }
}

void MdcTrackListBase::setD0Cut

(

double

d0Cut

)

[inline, inherited]

Definition at line 57 of file MdcTrackListBase.h.

References MdcTrackListBase::m_d0Cut.

{m_d0Cut = d0Cut;}//yzhang add 

void MdcTrackListBase::setPlot

(

int

plotFlag

)

[inline, inherited]

Definition at line 40 of file MdcTrackListBase.h.

References MdcTrackParams::lPlot, and MdcTrackListBase::tkParam.

{ tkParam.lPlot = plotFlag;};

void MdcTrackListBase::setPtCut

(

double

ptCut

)

[inline, inherited]

Definition at line 59 of file MdcTrackListBase.h.

References MdcTrackListBase::m_ptCut.

{m_ptCut = ptCut;}//yzhang add 2009-10-27

void MdcTrackListBase::setZ0Cut

(

double

z0Cut

)

[inline, inherited]

Definition at line 58 of file MdcTrackListBase.h.

References MdcTrackListBase::m_z0Cut.

{m_z0Cut = z0Cut;}//yzhang add 2010-05-21 

void MdcTrackListBase::store	(	RecMdcTrackCol *	,
		RecMdcHitCol *
	)			[inherited]

Definition at line 75 of file MdcTrackListBase.cxx.

References MdcTrackListBase::nTrack().

                                                                        {
  // ***********************************************************************
  int trackId = 0;
  for (int itrack = 0; itrack < nTrack(); itrack++) {
    MdcTrack* track = (*this)[itrack];
    //tkStat:  0,PatRec 1,MdcxReco 2,Tsf 3,CurlFinder
    int tkStat = 0;
    track->storeTrack(trackId, trackList, hitList, tkStat);
    ++trackId;
  }
  HepAListDeleteAll(*this);   // Discard the husks
  removeAll(); 
  return;
}

---

Member Data Documentation

double MdcTrackListBase::m_d0Cut = -999. [static, inherited]

Definition at line 60 of file MdcTrackListBase.h.

Referenced by createFromSegs(), MdcTrkRecon::initialize(), and MdcTrackListBase::setD0Cut().

int MdcTrackList::m_debug [private]

Definition at line 61 of file MdcTrackList.h.

Referenced by createFromSegs(), and finishHelix().

double MdcTrackListBase::m_ptCut = -999. [static, inherited]

Definition at line 62 of file MdcTrackListBase.h.

Referenced by createFromSegs(), MdcTrkRecon::initialize(), and MdcTrackListBase::setPtCut().

double MdcTrackListBase::m_z0Cut = -999. [static, inherited]

Definition at line 61 of file MdcTrackListBase.h.

Referenced by createFromSegs(), MdcTrkRecon::initialize(), and MdcTrackListBase::setZ0Cut().

MdcTrackParams MdcTrackListBase::tkParam [protected, inherited]

Definition at line 72 of file MdcTrackListBase.h.

Referenced by MdcTrackListBase::arbitrateHits(), MdcTrackListCsmc::createFromSegs(), createFromSegs(), MdcTrackListCsmc::finish3d(), finishCircle(), finishHelix(), MdcTrackListBase::MdcTrackListBase(), MdcTrackListBase::newParams(), pickHits(), and MdcTrackListBase::setPlot().

---