MdcTrackList Class Reference

#include <MdcTrackList.h>

Inheritance diagram for MdcTrackList:

List of all members.

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

---

Constructor & Destructor Documentation

MdcTrackList::MdcTrackList (  const MdcTrackParams &  tkPar  )

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

MdcTrackList::~MdcTrackList (   )

{}

MdcTrackList::MdcTrackList (  const MdcTrackList &   )  [private]

MdcTrackList::MdcTrackList (  const MdcTrackParams &  tkPar  )

MdcTrackList::~MdcTrackList (   )

MdcTrackList::MdcTrackList (  const MdcTrackList &   )  [private]

---

Member Function Documentation

int MdcTrackListBase::arbitrateHits (   )  [inherited]

int MdcTrackListBase::arbitrateHits (   )  [inherited]

{ //************************************************************************* // 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 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 if (8 == tkParam.lPrint) { std::cout<< "=remove above inactive "<<std::endl; } TrkFundHit* hit = const_cast<TrkFundHit*> (ihit->hit()); if(hitList->removeHit(hit)){ hit->setUnusedHit(&(*ihit)); } if(ihit == hitList->end()) break; --ihit;//be careful of the iterator, yzhang nDeleteInLayer[layer]++; } 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 for (int i=0;i<43;i++){ if (8 == tkParam.lPrint) { std::cout<<i<<" nHitInLayer "<<nHitInLayer[i] <<" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl; } if(nHitInLayer[i]>0 && (nHitInLayer[i]-nDeleteInLayer[i])==0) { nHitDeleted++; if (8 == tkParam.lPrint) { cout << "rec hits have been deleted in this layer"<<std::endl; } } } if (8 == tkParam.lPrint) { cout << " nHitDeleted = "<<nHitDeleted<<std::endl; } //yzhang add nHitDeleted cut 2010-09-13 // remove track if # not Active if (nHitDeleted >= tkParam.nHitDeleted) { if (8 == tkParam.lPrint) { cout << " nHitDeleted "<<nHitDeleted<<" > "<<tkParam.nHitDeleted <<" Killing tkNo " << itrack << endl; } 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()); if(hitList->removeHit(hit)){ alink->hit()->setUnusedHit(alink);//yzhang 2009-11-03 open } //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()); if(hitList->removeHit(hit)){ alink->hit()->setUnusedHit(alink);//yzhang 2009-11-03 open } //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.

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

Implements MdcTrackListBase. { // ************************************************************************* //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 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<<"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); #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); if(3==tkParam.lPrint){cout<<"Final fit2:"<<endl<<trialTrack->track()<< endl; } #ifdef MDCPATREC_TIMETEST TAU_PROFILE_STOP(timer7); #endif } // end if (success -- st) if(3==tkParam.lPrint){cout<<"Final fit3:"<<endl<<trialTrack->track()<< endl; } if (3 == m_debug){ dumpHelix(trialTrack); cout<<"finishHelix success? " << success<<"\n"; } if (!success) continue; if ( (m_z0Cut > Constants::epsilon) && (fabs(par.z0()) > m_z0Cut) ){ if (tkParam.lPrint>1)cout<<__FILE__ <<" Killing track by z0:" <<par.z0()<<">"<<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.

void MdcTrackList::dropMultiHotInLayer (  const MdcTrack *  tk  )

Reimplemented from MdcTrackListBase. { 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 *   )

void MdcTrackList::dumpAxCombine (  const MdcTrack *   )

{ if (trialTrack){ trialTrack->track().printAll(cout);//yzhang debug cout << "combined hits " << endl;//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; } }

void MdcTrackList::dumpAxFill (  const MdcTrack *   )

void MdcTrackList::dumpAxFill (  const MdcTrack *   )

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

void MdcTrackList::dumpCircle (  const MdcTrack *   )

void MdcTrackList::dumpCircle (  const MdcTrack *   )

{ if(trialTrack == NULL) return; if(trialTrack)trialTrack->track().printAll(cout); // TrkExchangePar par = trialTrack->track().fitResult()->helix(0.0); double omega,r,pt; if (!trialTrack->track().fitResult()) return; 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 //m_nHitCir = trialTrack->track().hits()->hotList().nHit(); std::cout << "hitList after finishCircle, nHit=" << trialTrack->track().hits()->hotList().nHit()<< std::endl;//yzhang debug TrkHotList::hot_iterator hotIter= trialTrack->track().hits()->hotList().begin(); while (hotIter!=trialTrack->track().hits()->hotList().end()) { //m_nHitCir = ((MdcHit*)(hotIter->hit()))->layernumber();//yzhang debug cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber() <<","<<((MdcHit*)(hotIter->hit()))->wirenumber() <<":"<<hotIter->isActive()<<") "; //double tphi = gm->Wire(((MdcHit*)(hotIter->hit()))->layernumber(), // ((MdcHit*)(hotIter->hit()))->wirenumber())->phiE(); hotIter++; //std::cout <<"in MdcTrackList Circle tphiEnd "<< tphi << std::endl;//yzhang debug //m_lwPhi = tphi; } cout <<endl; }

void MdcTrackList::dumpD0 (  const TrkExchangePar &   )

void MdcTrackList::dumpD0 (  const TrkExchangePar &   )

{ //yzhang hist //m_hd0->fill(par.d0()); //m_d0 = par.d0(); // m_tuple1->write();//yzhang hist std::cout << " MdcTrackList::par.d0() " << par.d0()<<"\n";//yzhang debug //zhangy hist }

void MdcTrackList::dumpHelix (  const MdcTrack *   )

void MdcTrackList::dumpHelix (  const MdcTrack *   )

{ std::cout << "hitList after finishHelix, nHit=" << trialTrack->track().hits()->hotList().nHit()<< std::endl;//yzhang debug 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; }

void MdcTrackList::dumpSeed (  const MdcSeg *  seed, bool  goodOnly )

void MdcTrackList::dumpSeed (  const MdcSeg *  seed, bool  goodOnly )

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

void MdcTrackList::dumpStCombine (  const MdcTrack *   )

void MdcTrackList::dumpStCombine (  const MdcTrack *   )

{ std::cout << "hitList after StCombine" << std::endl;//yzhang debug 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() <<":"<<hotIter->isActive() <<") "; hotIter++; tmplay=layer; } cout <<endl; }

void MdcTrackList::dumpStFill (   )

void MdcTrackList::dumpStFill (   )

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

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

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

{ //************************************************************************ TrkRecoTrk& trk = mdcTrk.track(); if(9==tkParam.lPrint){ std::cout << __FILE__ << " 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, ?? int nHits = pickHits(&mdcTrk, map, gm, lcurler); if (nHits <= 0) { return success; }//yzhang FIXME 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 *  )

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

{ //*********************************************************************** 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, ?? int nHits = pickHits(&mdcTrk, map, gm, lcurler);//yzhang add 2010-05-10 //if (nHits <= 0) { return success; }//yzhang FIXME 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 if(3==tkParam.lPrint){ cout << "Final fit1: " << endl << trk << endl; } return success; }

void MdcTrackListBase::newParams (  const MdcTrackParams &  tkPar  )  [inherited]

void MdcTrackListBase::newParams (  const MdcTrackParams &  tkPar  )  [inherited]

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

int MdcTrackListBase::nTrack (   )  const [inline, inherited]

{return length();}

int MdcTrackListBase::nTrack (   )  const [inline, inherited]

{return length();}

MdcTrackList& MdcTrackList::operator= (  const MdcTrackList &   )  [private]

MdcTrackList& MdcTrackList::operator= (  const MdcTrackList &   )  [private]

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

int MdcTrackList::pickHits (  MdcTrack *  trk, const MdcHitMap *  map, const MdcDetector *  gm, bool  pickAm = true )

!! change to using omega itself { /***************************************************************************/ // 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 if(6==tkParam.lPrint){ cout << "Before pickHits: " << endl << trk << 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 // *** 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; //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 delete //**** 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" << std::endl; } 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" << std::endl; } //stereo layer pick hit phi wider double maxInactiveResid = tkParam.maxInactiveResid; //std::cout<< __FILE__ <<" maxInactiveResid "<< maxInactiveResid <<" pickHitPhiFactor "<<tkParam.pickHitPhiFactor<< std::endl; if ( 0 != layer->view()) maxInactiveResid *= (tkParam.pickHitPhiFactor*fabs(layer->dPhiz()));//FIXME yzhang 2010-05-18 , original 5 //std::cout<< " facotr "<<tkParam.pickHitPhiFactor<< " dPhiz "<<layer->dPhiz() //<< " factor*dPhiz "<<layer->dPhiz()*tkParam.pickHitPhiFactor<< std::endl; if (curv >= 0.0) { 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; } 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<<__FILE__ << " phiLow "<<phiLow << " phiHigh "<<phiHigh << std::endl; std::cout<<__FILE__ << " wireLow "<<wireLow << " wireHigh "<<wireHigh << std::endl; } // *** Loop through the predicted hit wires for (int thisWire = wireLow; thisWire <= wireHigh; thisWire++) { int corrWire = mdcWrapWire(thisWire, nCell); thisHit = map->hitWire(layer->layNum(), corrWire); double tof = 0.; double z = 0.; double driftDist = 0.; if(6==tkParam.lPrint){ if(thisHit != 0) thisHit->print(std::cout); else std::cout << " test "<<layer->layNum()<<" "<<corrWire<< std::endl; } // Calculate expected drift distance double delx, dely; double resid = 0., predDrift = 0.; int ambig = 0; const MdcHitOnTrack *alink = 0; 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); } else { // look for existing hit //if (thisHit->nUsedHits()==0) break;//yzhang temp //yzhang CHANGE if(thisHit->getHitOnTrack(&(trk->track())) ==0){ if(6==tkParam.lPrint)std::cout << "getHit==0 " << std::endl;//yzhang debug alink = 0;//yzhang temp }else{ alink = thisHit->getHitOnTrack(&(trk->track()))->mdcHitOnTrack();//yzhang temp } if(6==tkParam.lPrint) { cout << " Fetched existing hot " << endl;//yzhang debug } if (alink == 0 || pickAm) { if ((!tkStatus->is2d()) && layer->view() != 0){ double rc = 9999.; if(fabs(par.omega())>Constants::epsilon) rc = 1./par.omega(); //std::cout<< par.omega()<<" "<< rc << std::endl; double rw = layer->rMid(); double alpha = acos(1 - rw*rw/(2*rc*rc)); //std::cout<< rw <<" "<< 1 - rw*rw/(2*rc*rc)<<" alpha " <<alpha << std::endl; double fltLen = rw; if(fabs(1 - rw*rw/(2*rc*rc))<1) fltLen = rc * alpha; z = par.z0() + fltLen* par.tanDip(); //std::cout<<"z "<<z<<" fltLen "<<fltLen<< " tanl "<<par.tanDip() << std::endl; 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; //std::cout<< "new x "<<x <<" y: "<< y << std::endl; //std::cout<< layer->xWire(corrWire) <<" old y: "<< layer->yWire(corrWire) << std::endl; }else{ delx = -d0 * sinPhi0 - thisHit->x(); dely = d0 * cosPhi0 - thisHit->y(); } predDrift = cosPhi0 * dely - sinPhi0 * delx + curv * (delx * delx + dely * dely); if(6==tkParam.lPrint){ if (alink == 0) cout << "hit not on track: pred drift dist = " << predDrift << endl;//yzhang debug else cout << "picking all " << predDrift << endl;//yzhang debug } // predDrift = predDrift * (1. - curv() * predDrift); ambig = (predDrift >= 0) ? 1 : -1; if (hasCurled) ambig = -ambig; //std::cout<<bunchTime<<" bunchTime "<< std::endl; double entranceAngle=0.; driftDist = thisHit->driftDist(tof,ambig,entranceAngle,0.,z); //FIXME //std::cout<< thisHit->driftDist(bunchTime,ambig) <<" new drift "<< driftDist << std::endl; if (alink == 0) { // FIXME: is this ambig here VVVVV OK for incoming tracks? //resid = ambig * fabs(thisHit->driftDist(bunchTime,ambig)) - predDrift; resid = ambig * fabs(driftDist) - predDrift; aresid = fabs(resid); if (aresid > tkParam.maxInactiveResid ) ambig = 0; } if(6==tkParam.lPrint) cout << " pickAm ambig " << ambig << " drift "<<driftDist<< " predDrift "<<predDrift<<" resid "<<resid<< endl;//yzhang debug } else { ambig = alink->ambig(); if(6==tkParam.lPrint) cout << " alink==0 ambig " << ambig<< " drift "<<driftDist<< " predDrift "<<predDrift<<" resid "<<resid<< endl;//yzhang debug } } // Cases // (1) No good hit, but track near cell-edge => do nothing and continue if(6==tkParam.lPrint) std::cout << __FILE__ << " predDrift "<<fabs(predDrift)<<">? goodDriftCut "<<3. * goodDriftCut<< std::endl; if (ambig == 0 && fabs(predDrift) > 3.*goodDriftCut){//yzhang 2009-11-03 add 3factor if(6==tkParam.lPrint) cout << " (1) No good hit, but track near cell-edge " << endl;//yzhang debug 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; double residCut = tkParam.maxActiveSigma * factor * sigma; if(6==tkParam.lPrint){ std::cout<< "aresid "<<aresid<<"<?" <<residCut <<" maxActiveSigma "<<tkParam.maxActiveSigma <<" factor="<<factor <<" sigma "<<sigma<< std::endl; } if(fabs(par.tanDip())<2) factor = (2-par.tanDip())/2 * factor; if (alink == 0 && (aresid <= residCut ) ) { //if (alink == 0) int isActive = 1; if (aresid > residCut) isActive = 0; if(6==tkParam.lPrint) { if (aresid > residCut) std::cout << "notACT, resid "<<aresid<<" > " << 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); //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<< "appendHot " << std::endl; alink = (MdcHitOnTrack*) hitList->appendHot(&tempHot); } if (!localHistory.member(const_cast<MdcHitOnTrack*>(alink))) { localHistory.append(const_cast<MdcHitOnTrack*>(alink)); if (hasCurled) trk->setHasCurled(); nFound++; 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; } if(6==tkParam.lPrint) cout << " (2) Hit found " << endl;//yzhang debug } // (3) No hit found; if beyond known good region, should we continue? else if (ambig == 0 && reachedLastInput) { 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) { //std::cout<< __FILE__ <<" "<< __LINE__ <<" break by lContinue "<< 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 (!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; } } } //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 if(6==tkParam.lPrint){ cout << " localHistory= 0 " << endl; } if(6==tkParam.lPrint){ cout << "After pickHits found " << nFound <<" hits"<< endl; hitList->hotList().print(std::cout); } return nFound; }

void MdcTrackListBase::plot (   )  const [inherited]

void MdcTrackListBase::plot (   )  const [inherited]

{ //************************************************************************* 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, and MdcTrackListCsmc.

void MdcTrackListBase::remove (  MdcTrack *  atrack  )  [inherited]

Reimplemented in MdcTrackListCsmc, and MdcTrackListCsmc. { //-------------------------------------------------------------------- if (atrack != 0) { HepAList<MdcTrack>::remove( atrack ); delete atrack; } }

void MdcTrackListBase::setD0Cut (  double  d0Cut  )  [inline, inherited]

{m_d0Cut = d0Cut;}//yzhang add

void MdcTrackListBase::setD0Cut (  double  d0Cut  )  [inline, inherited]

{m_d0Cut = d0Cut;}//yzhang add

void MdcTrackListBase::setPlot (  int  plotFlag  )  [inline, inherited]

{ tkParam.lPlot = plotFlag;};

void MdcTrackListBase::setPlot (  int  plotFlag  )  [inline, inherited]

{ tkParam.lPlot = plotFlag;};

void MdcTrackListBase::setPtCut (  double  ptCut  )  [inline, inherited]

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

void MdcTrackListBase::setPtCut (  double  ptCut  )  [inline, inherited]

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

void MdcTrackListBase::setZ0Cut (  double  z0Cut  )  [inline, inherited]

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

void MdcTrackListBase::setZ0Cut (  double  z0Cut  )  [inline, inherited]

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

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

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

{ // *********************************************************************** 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]

int MdcTrackList::m_debug [private]

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

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

MdcTrackParams MdcTrackListBase::tkParam [protected, inherited]

---

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

• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcTrackList.h
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Reconstruction/MdcPatRec/MdcTrkRecon/MdcTrkRecon-00-03-31/MdcTrkRecon/MdcTrackList.h
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Reconstruction/MdcPatRec/MdcTrkRecon/MdcTrkRecon-00-03-31/src/MdcTrackList.cxx

---