AlgTrackSRList Class Reference

#include <AlgTrackSRList.h>

Inheritance diagram for AlgTrackSRList:
AlgBase

List of all members.

Public Member Functions

 AlgTrackSRList ()
virtual ~AlgTrackSRList ()
virtual void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void Trace (const char *c) const

Private Member Functions

void AddBestPlaneClusterToTrack (TrackClusterSR *besttc, Track2DSR *track)
Int_t AddClustersToTracks (Track2DSRItr &trackItr, TrackClusterSRItr &planeClusterItr, Int_t direction, Int_t iterate, std::vector< Int_t > &hitPlanes, Int_t trk2dmin)
Bool_t CheckForBadClusters (TrackClusterSR *besttc, Track2DSR *track, Int_t &planesBefore, Int_t direction, Int_t trk2dmin)
void CleanUp (TObjArray *trackList, TObjArray *clusterList)
void DeleteTwinTracks (CandTrackSRHandleItr &candTracks1, CandTrackSRHandleItr &candTracks2, TObjArray *candTracks, Int_t uMaxTrack, Int_t vMaxTrack, Int_t nmax3dtrk, CandHandle &ch)
void FillInGaps (Track2DSRItr &trackItr, TrackClusterSRItr &clusterItr, Int_t direction)
Int_t FillTrackList (TrackClusterSRItr &clusterItr, TObjArray *trackList, TObjArray *seedClusters, Double_t houghSlope, Double_t houghIntercept, Int_t direction, Int_t iterate, Int_t &begPlane)
void Find2DTrackEndPoints (CandStripHandleItr &stripItr, Int_t &vtxPlane, Int_t &endPlane, Int_t &numPlanes)
Bool_t PlaneIsActive (Int_t plane, Float_t u, Float_t v, Float_t tol)
Int_t FindNumSkippedPlanes (Int_t currentPlane, Int_t prevPlane, std::vector< Int_t > &hitPlanes, Int_t direction)
Int_t FindNumSkippedPlanes (Int_t currentPlane, Int_t prevPlane, Float_t prevTPos, Float_t prevZPos, Double_t trackSlope, Int_t direction, PlaneView::PlaneView_t view)
Int_t FindNumSkippedPlanes (Int_t currentPlane, Int_t prevPlane, Float_t PrevUPos, Float_t PrevVPos, Float_t prevZPos, Double_t trackSlopeU, Double_t trackSlopeV, Int_t direction, PlaneView::PlaneView_t view)
Int_t FindTimingDirection (CandStripHandleItr &allStripItr, Float_t uTrackSlope, Float_t uTrackIntercept, Float_t vTrackSlope, Float_t vTrackIntercept)
Int_t FindTimingDirection (TrackClusterSRItr &clusterItr)
void FormCandTrackSR (Track2DSRItr &uTrackItr, Track2DSRItr &vTrackItr, TObjArray *candTracks, CandHandle &ch, CandSliceHandle *slice, AlgHandle ah, CandContext cxx, Int_t direction)
Int_t IdentifyBadTracks (Track2DSRItr &trackItr, Int_t iterate, Int_t trk2dmin)
Bool_t IsBestClusterInPlane (TrackClusterSR *planeCluster, TrackClusterSR *prevCluster, Int_t numSkippedHit, Int_t nHit, Int_t planesBefore, Double_t trackSlope, Track2DSR *track, Double_t *residParams, Int_t dplane, Bool_t isSpectrometerPlane)
Int_t LookForBadClusters (TrackClusterSRItr &clusterItr, TObjArray *trackList, Int_t direction)
void MakeTrackClusters (TObjArray *trackClusterList, CandStripHandleItr &stripItr, CandTrackSRListHandle &cth, Int_t expectedStrips, Int_t direction)
void MarkUsedClusters (Track2DSRItr &trackItr, TrackClusterSRItr &clusterItr, Int_t iterate, Int_t nmaxtrack, Double_t houghSlope, Double_t houghIntercept)
Int_t RemoveTrackSubsets (Track2DSRItr &trackItr1, Track2DSRItr &trackItr2, TObjArray *trackList)
void SpectrometerTracking (TObjArray *trackList, CandTrackSRListHandle &cth, CandSliceHandle *slice, CandSliceHandle *dupSlice, CandContext cx)
void MergeTracks (TObjArray *tracklist, Int_t direction)
void RemoveUnusedSpectStrips (CandTrackSRListHandle &cth, CandStripListHandle *striplist)
void RemoveStripsInSlice (CandTrackSRListHandle &cth, CandSliceListHandle *slicelist)

Private Attributes

Detector::Detector_t fDetector
SimFlag::SimFlag_t fSimFlag
Int_t fMapIsWide [1000]
Int_t fSliceVtxUPlane
Int_t fSliceVtxVPlane
Int_t fSliceEndUPlane
Int_t fSliceEndVPlane
Int_t fSliceVtxPlane
Int_t fSliceEndPlane
Int_t fSliceNUPlanes
Int_t fSliceNVPlanes
Int_t fSliceNPlanes
Double_t fSliceDzDs
VldContext fVldc
PlaneOutline fPL
Int_t fParmMaxNStrip
Int_t fParmMaxNExtraStrip
Int_t fParmTrk2DPlnEnd
Double_t fParmTrk2DWin0
Double_t fParmHitNTime
Double_t fParmHitTime
Int_t fParmTrk2DNSkip
Int_t fParmTrk2DNSkipHit
Double_t fParmTrk2DAlpha
Double_t fParmTrk2DMaxResid
Int_t fParmTrk2DNSkipRemove
Int_t fParmTrkNPlaneGoodDir
Double_t fParmTrk2DDirCosNSigma
Double_t fParmTrk2DDirCosError2
Int_t fParmTrk2DEndPlaneDiff
Int_t fParmMinSingleHit
Int_t fParmSingleHitDef
Int_t fParmTrk2DNHough0
Double_t fParmTrk2DLinA0
Double_t fParmTrk2DLinB0
Int_t fParmTrk2DNHough
Double_t fParmTrk2DLinA
Double_t fParmTrk2DLinB
Double_t fParmTrk2DNSigmaA
Int_t fParmTrk2DNSeed
Int_t fParmTrk2DNContiguous
Double_t fParmTrk2DHitFraction
Double_t fParmTrk2DTwinMatchFraction
Int_t fParmTrk2DSubsetNHit
Int_t fParmTrk2DSubsetDHit
Int_t fParmDiffViewBegPlaneMatch
Int_t fParmDiffViewEndPlaneMatch
Double_t fParmDiffViewPlaneOverlap
Double_t fParmDiffViewTimeMatch
Double_t fParmTrk3DTwinFrac
Int_t fParmTrkClsNSkip
Int_t fParmTrk3DTrackMax
Int_t fParmIsCosmic
Double_t fParmMinStripPulseHeight
Double_t fParmMinClusterPulseHeight
Int_t fParmUseWideClusters
Int_t fParmUseShowerLikePlanes
Double_t fParmMaxTimingResid
Double_t fParmMisalignmentError
Int_t fParmHoughMinBin
Double_t fParmHoughMinFrac
Double_t fParmHoughMinFracZoom
Double_t fParmHoughMinInterBinSize
Double_t fParmExtrapError
CandDigitListHandlefCDLHnew
CandStripListHandlefCSLHnew
CandSliceListHandlefCSlcLHnew
Float_t fHoughUSlope
Float_t fHoughVSlope
Float_t fHoughUIntercept
Float_t fHoughVIntercept

Detailed Description

Definition at line 40 of file AlgTrackSRList.h.


Constructor & Destructor Documentation

AlgTrackSRList::AlgTrackSRList (  ) 

Definition at line 176 of file AlgTrackSRList.cxx.

00177 {
00178   for(Int_t i = 0; i<1000; ++i){
00179     fMapIsWide[i] = 0;
00180   }
00181 }

AlgTrackSRList::~AlgTrackSRList (  )  [virtual]

Definition at line 184 of file AlgTrackSRList.cxx.

00185 {
00186 }


Member Function Documentation

void AlgTrackSRList::AddBestPlaneClusterToTrack ( TrackClusterSR besttc,
Track2DSR track 
) [private]

Definition at line 2888 of file AlgTrackSRList.cxx.

References Track2DSR::Add(), fParmTrk2DAlpha, Track2DSR::GetCluster(), Track2DSR::GetForwardSlope(), Track2DSR::GetLast(), TrackClusterSR::GetTPos(), and TrackClusterSR::GetZPos().

Referenced by AddClustersToTracks().

02890 {
02891   
02892   Int_t iclosest=-1;
02893   Float_t dzclosest=1e6;
02894   Int_t ilast = track->GetLast();
02895   TrackClusterSR *tc0 = track->GetCluster(ilast);
02896   for (Int_t i=0; i<track->GetLast(); i++){
02897     tc0 = track->GetCluster(i);
02898     if(TMath::Abs(tc0->GetZPos()-besttc->GetZPos())<dzclosest){
02899       iclosest=i;
02900       dzclosest=TMath::Abs(tc0->GetZPos()-besttc->GetZPos());
02901     }
02902   }
02903   Double_t slope0=0;
02904   if(iclosest>=0){
02905     slope0 = track->GetForwardSlope(iclosest);
02906   }
02907   Double_t slope = (besttc->GetTPos()-tc0->GetTPos()); 
02908   slope /= (besttc->GetZPos()-tc0->GetZPos());
02909   
02910   //weight the slope of the this cluster using the previous clusters in the track
02911   if (slope0!=0) {
02912     slope *= (1.-fParmTrk2DAlpha);
02913     slope += fParmTrk2DAlpha*slope0;
02914   }
02915   track->Add(besttc,slope);
02916   return;
02917 }

Int_t AlgTrackSRList::AddClustersToTracks ( Track2DSRItr &  trackItr,
TrackClusterSRItr &  planeClusterItr,
Int_t  direction,
Int_t  iterate,
std::vector< Int_t > &  hitPlanes,
Int_t  trk2dmin 
) [private]

Definition at line 2152 of file AlgTrackSRList.cxx.

References AddBestPlaneClusterToTrack(), CheckForBadClusters(), Track2DSR::Compress(), fDetector, fHoughUIntercept, fHoughUSlope, fHoughVIntercept, fHoughVSlope, FindNumSkippedPlanes(), Track2DSR::fIterate, fParmTrk2DNContiguous, fParmTrk2DNSkip, fParmTrk2DNSkipHit, fParmTrk2DNSkipRemove, Track2DSR::GetBegPlane(), Track2DSR::GetCluster(), Track2DSR::GetForwardSlope(), Track2DSR::GetLast(), TrackClusterSR::GetPlane(), VHS::GetPlane(), TrackClusterSR::GetPlaneView(), TrackClusterSR::GetTPos(), TrackClusterSR::GetZPos(), Track2DSR::IsBad(), IsBestClusterInPlane(), TrackClusterSR::IsValid(), Msg::kDebug, Detector::kNear, PlaneView::kU, Msg::kVerbose, MSG, and Track2DSR::RemoveAt().

Referenced by RunAlg().

02156 {
02157 
02158   trackItr.ResetFirst();
02159 
02160   Int_t currentPlane = 0;
02161   Int_t hitPlanesSkippedInTrack = 0;
02162   TrackClusterSR *besttc=0;  //place holder of the cluster to add
02163   TrackClusterSR *prevCluster=0;  //previous cluster in track to current plane
02164   TrackClusterSR *planeCluster=0; //cluster of current plane we are looking at
02165   Track2DSR *track = 0; //current track
02166 
02167   //set the direction of your iterators
02168   if(direction == 1) planeClusterItr.ResetFirst();
02169   else if(direction == -1) planeClusterItr.ResetLast();
02170  
02171   Int_t planesBefore = 0;
02172 
02173   //numSkippedActive = # of active planes we cant ignore between
02174   //the current plane and the plane of the last cluster in the current 
02175   //track.  numSkippedHit = # of hit planes in this event skipped
02176   Int_t numSkippedActive = 0;
02177   Int_t numSkippedHit = 0;
02178 
02179   //this is the slope of the track at the cluster you are looking at
02180   Double_t trackSlope = 0.;
02181 
02182   //hold the values of the fit for the clusters in the following array
02183   //residParams[0] is the best residual without using signal weights, 
02184   //residParams[1] is the best residual using signal weights
02185   Double_t residParams[] = {1.e6,1.e6};
02186 
02187   //number of planes between current cluster and previous one already in track
02188   Int_t deltaPlane = 0;
02189   Int_t trackCtr = 0;
02190   
02191   //# of planes between the one with the last cluster added to the track
02192   //and the plane in the track holding prevCluster
02193   Int_t nHit = 0;
02194 
02195   //make an array to hold the index of hitPlanes for a given plane
02196   Int_t planeIndex[500];
02197   for(Int_t i = 0; i < 500; ++i){
02198     planeIndex[i] = 0;
02199   }
02200   for(Int_t ii = 0; ii < static_cast<Int_t>(hitPlanes.size()); ++ii){
02201    planeIndex[hitPlanes[ii]] = ii;
02202   
02203   }
02204   Int_t index = 0;
02205 
02206   //loop over the tracks in the list - this way you add clusters from
02207   //the current plane to each track and save cpu time vs looping over
02208   //each track and then every plane in the list
02209   trackCtr = 0;
02210   while( (track = trackItr()) ){   
02211     MSG("TrackSR", Msg::kDebug) << "iterate = " << iterate 
02212                                 << " track->fIterate = " << track->fIterate 
02213                                 << " track # " << trackCtr 
02214                                 << " track beg plane = " << track->GetBegPlane() 
02215                                 << " track is bad = " << (Int_t)track->IsBad() 
02216                                 << endl;
02217     index = planeIndex[track->GetBegPlane()];  
02218      
02219     //if this track was created in this iteration and was not
02220     //previously marked as a bad track
02221     if( track->fIterate == iterate && !track->IsBad() ){
02222       
02223       MSG("AlgTrackSRList", Msg::kVerbose) << "on a good track in iteration " 
02224                                          << iterate << endl;
02225       planesBefore = 0;
02226       planeClusterItr.ResetFirst();
02227 
02228       //loop over the planes hit in this slice
02229       for(Int_t plnIndex=index; plnIndex<static_cast<Int_t>(hitPlanes.size()); 
02230           ++plnIndex){
02231         
02232         currentPlane = hitPlanes[plnIndex];
02233         deltaPlane = 0;
02234         MSG("TrackSR", Msg::kDebug) << "currentPlane = " << currentPlane << endl;                       
02235         numSkippedActive = 0;
02236         numSkippedHit = 0;
02237         
02238         //have to reset planesBefore for each new track
02239         planesBefore = 0;
02240                 
02241         //reset besttc because you are on a new plane
02242         besttc = 0;
02243         
02244         // make sure that the current plane is downstream of 
02245         // the first cluster in the list of clusters for this time through the finder
02246 
02247         if(currentPlane*direction
02248            >(direction*track->GetCluster(track->GetLast())->GetPlane())){ 
02249               
02250           //loop from the current last cluster in the track to the beginning of the
02251           //track if necessary to see if the cluster of the next plane matches
02252           //with any of them
02253 
02254           MSG("TrackSR", Msg::kDebug) << "planesBefore = " 
02255                                       << planesBefore << endl;
02256 
02257           //loop backwards over the clusters in the track
02258           while( planesBefore<= fParmTrk2DNSkipRemove 
02259                  && track->GetLast()-planesBefore>=0
02260                  && !besttc){
02261 
02262             //get the cluster you are interested in
02263             prevCluster = track->GetCluster(track->GetLast()-planesBefore);
02264             
02265             MSG("TrackSR", Msg::kDebug) << track->GetForwardSlope(track->GetLast()-planesBefore)
02266                                         << " " << planesBefore << " " 
02267                                         << track->GetLast()+1 << endl;
02268 
02269             trackSlope = track->GetForwardSlope(track->GetLast()-planesBefore);
02270             Float_t upos=0;
02271             Float_t vpos=0;
02272             Float_t uslope=0;
02273             Float_t vslope=0;
02274             if(prevCluster->GetPlaneView()==PlaneView::kU){
02275               upos = prevCluster->GetTPos();
02276               vpos = fHoughVIntercept + prevCluster->GetZPos()*fHoughVSlope;
02277               uslope=track->GetForwardSlope(track->GetLast()-planesBefore);
02278               vslope=fHoughVSlope;
02279             }
02280             else{
02281               vpos = prevCluster->GetTPos();
02282               upos = fHoughUIntercept + prevCluster->GetZPos()*fHoughUSlope;
02283               vslope=track->GetForwardSlope(track->GetLast()-planesBefore);
02284               uslope=fHoughUSlope;
02285             }
02286             numSkippedActive=FindNumSkippedPlanes(currentPlane,
02287                                                   prevCluster->GetPlane(),
02288                                                   upos,vpos,
02289                                                   prevCluster->GetZPos(),
02290                                                   uslope,
02291                                                   vslope,
02292                                                   direction,
02293                                                   prevCluster->GetPlaneView());
02294             
02295             deltaPlane = TMath::Abs(prevCluster->GetPlane()-currentPlane);
02296             
02297             MSG("TrackSR", Msg::kDebug) <<"current plane " 
02298                                         << currentPlane 
02299                                         <<" prev plane " 
02300                                         << prevCluster->GetPlane() 
02301                                         <<" active skipped = " 
02302                                         << numSkippedActive
02303                                         <<" deltaPlane = " 
02304                                         << deltaPlane << endl;
02305             
02306             //the +1 below accounts for the fact that get last
02307             //returns a number starting the count in the array at 0
02308             nHit = track->GetLast()-planesBefore+1;
02309             
02310             //if the distance between the plane of the previous cluster and
02311             //the current plane seems reasonable, then look for clusters 
02312             //to add from the current plane.  make sure you arent looking 
02313             //at clusters from the same plane though
02314             if(numSkippedActive<=fParmTrk2DNSkip
02315                && (nHit>=fParmTrk2DNContiguous || numSkippedActive==0)
02316                && currentPlane != prevCluster->GetPlane() ){
02317               
02318               //reset the residParams because you are on a new plane
02319               residParams[0] = 1.e6;
02320               residParams[1] = 1.e6;
02321 
02322              
02323               //slice the planeClusterItr to get only those clusters on the current plane
02324               planeClusterItr.GetSet()->Slice();
02325               planeClusterItr.GetSet()->Slice(currentPlane);     
02326               //loop over the clusters in the current plane
02327               while( planeClusterItr.IsValid() ){
02328                 planeCluster = planeClusterItr.Ptr();
02329                 MSG("TrackSR", Msg::kDebug) << " prev cluster plane " << prevCluster->GetPlane() << " tpos  " << prevCluster->GetTPos() << " valid  " << prevCluster->IsValid()  << " cluster " << planeCluster->GetPlane() << " tpos  " << planeCluster->GetTPos() << " valid " << planeCluster->IsValid() << endl;
02330                                 
02331                 Bool_t isSpectrometerPlane = false; 
02332                 if(fDetector==Detector::kNear 
02333                    && (planeCluster->GetPlane()>=121)){
02334                   isSpectrometerPlane=true;
02335                 }
02336                 if(planeCluster->IsValid()){
02337                   if(IsBestClusterInPlane(planeCluster,prevCluster, 
02338                                           numSkippedHit,
02339                                           nHit, planesBefore,trackSlope, track, residParams, 
02340                                           deltaPlane,isSpectrometerPlane)) besttc = planeClusterItr.Ptr();
02341                 }
02342                 planeClusterItr.Next();
02343               }//end loop over the clusters in plane currentPlane
02344               
02345               //check to see if there are any clusters which need to be removed
02346               //from the fit between the plane of planes before and that of 
02347               //planeCluster with the addition of planeCluster the 
02348               //CheckForBadClusters returns a true if it is ok to remove bad clusters
02349               //and zero's besttc if there were not - so if it is true,
02350               //then we add besttc to the track
02351               if(CheckForBadClusters(besttc, track, planesBefore, direction, 
02352                                      trk2dmin) ){
02353                 
02354                 MSG("TrackSR", Msg::kDebug) << " adding tc "
02355                                             << besttc->GetPlane()
02356                                             << " " 
02357                                             << besttc->GetTPos()
02358                                             << endl;
02359                 AddBestPlaneClusterToTrack(besttc,track);
02360                 
02361                 //reset hitPlanesSkippedInTrack because you just added a hit plane to the track
02362                 hitPlanesSkippedInTrack = 0;
02363               }
02364               //no besttc on this plane
02365               else if(direction*currentPlane
02366                       >direction*track->GetCluster(track->GetLast())->GetPlane())
02367                 ++hitPlanesSkippedInTrack;
02368               
02369             }//end if distances between planes is reasonable
02370               
02371             //the distances were not reasonable, so why waste any more 
02372             //time on this track by making the distance even larger?  
02373             //set planesBefore to 999999 to break out of the loop for 
02374             //this track
02375             else planesBefore = 999999;
02376             
02377             MSG("TrackSR", Msg::kDebug) << "planes before = " << planesBefore 
02378                                         << " track->GetLast() = " 
02379                                         << track->GetLast()+1
02380                                         << endl;
02381             ++planesBefore;
02382           }//end loop to look for a cluster on this plane
02383         }//end if currentPlane is downstream of the beginning plane for this loop  
02384         
02385         //if you have a stupid number for planes before, stop the for loop and move on to the next track
02386         if(planesBefore >= 999999) break;
02387       }//end loop over planes hit in slice
02388     }//end if track made in this iteration
02389 
02390     ++trackCtr;
02391     
02392   }//end loop over tracks
02393 
02394   trackItr.ResetFirst();
02395   // check each track for outlier hit at two track ends, and remove if found
02396   while( (track = trackItr()) ){
02397     if(track->GetLast()>2){
02398       Int_t gap1=TMath::Abs(track->GetCluster(0)->GetPlane()-track->GetCluster(1)->GetPlane())/TMath::Abs(track->GetCluster(2)->GetPlane()-track->GetCluster(1)->GetPlane());
02399       Int_t gap2=TMath::Abs(track->GetCluster(track->GetLast())->GetPlane()-track->GetCluster(track->GetLast()-1)->GetPlane())/TMath::Abs(track->GetCluster(track->GetLast()-1)->GetPlane()-track->GetCluster(track->GetLast()-2)->GetPlane());
02400       Int_t gapfactor=2;
02401       if(gap1>fParmTrk2DNSkipHit*gapfactor)track->RemoveAt(0);
02402       if(gap2>fParmTrk2DNSkipHit*gapfactor)track->RemoveAt(track->GetLast());
02403       track->Compress();
02404     }   
02405   }
02406   trackItr.ResetFirst();
02407   return hitPlanesSkippedInTrack;
02408 }

Bool_t AlgTrackSRList::CheckForBadClusters ( TrackClusterSR besttc,
Track2DSR track,
Int_t &  planesBefore,
Int_t  direction,
Int_t  trk2dmin 
) [private]

Definition at line 2802 of file AlgTrackSRList.cxx.

References Track2DSR::Compress(), fParmTrk2DHitFraction, fParmTrk2DMaxResid, fParmTrk2DNHough, Track2DSR::GetBegPlane(), Track2DSR::GetCluster(), Track2DSR::GetLast(), TrackClusterSR::GetPlane(), VHS::GetPlane(), TrackClusterSR::GetTPos(), TrackClusterSR::GetZPos(), Msg::kDebug, max, MSG, Track2DSR::RemoveAt(), and LinearFit::Unweighted().

Referenced by AddClustersToTracks().

02806 {
02807   
02808   Bool_t removetc(kTRUE);
02809 
02810   if( !besttc) return false;
02811   if(planesBefore<=0 ) return true;
02812   
02813    MSG("TrackSR", Msg::kDebug) << " ----------checking besttc on plane------------ " << besttc->GetPlane()
02814                               << endl;
02815 
02816    //get the number of clusters (ie hit planes) in the track.
02817    //add one to track->GetLast() because that returns a number starting
02818    //from 0 - its an array count after all, and add one to planesBefore
02819    //to account for besttc just found but not yet added
02820    Int_t nhit = track->GetLast()+1-planesBefore+1;
02821    
02822    //assume alternate views when finding dplane - the number of planes
02823    //between the first plane in the track and the plane of the current cluster.
02824    //not the best way since it has problems in tracks that jump the super module
02825    //gap, but good enough for now
02826    Int_t dplane = direction*(besttc->GetPlane()-track->GetBegPlane())/2+1;
02827    
02828    MSG("TrackSR", Msg::kDebug) << "nhit = " << nhit << " dplane = " << dplane
02829                                << "/" << fParmTrk2DHitFraction << " trk2dmin = "
02830                                << trk2dmin << " 2DNHough = " << fParmTrk2DNHough
02831                                << endl;
02832    
02833    //if the # of hits in the track is less than the hit fraction * dplane
02834    //you cant remove any of the hits - you need them all
02835    //if the # of hits is less than the minimum you cant get rid of any 
02836    if ((1.*nhit)<fParmTrk2DHitFraction*dplane || nhit<trk2dmin) {
02837      removetc = kFALSE;
02838      planesBefore = 9999999;
02839    }
02840    
02841    //if you only have enough hits to do hough tracks, see how well
02842    //the current hit fits in with the rest
02843    else if(nhit<=fParmTrk2DNHough){
02844      Double_t *xfit = new Double_t[nhit];
02845      Double_t *yfit = new Double_t[nhit];
02846      Double_t parm[2];
02847      TrackClusterSR *tctmp = 0;
02848      for (Int_t ifit=0; ifit<nhit; ++ifit){
02849        tctmp = track->GetCluster(ifit);
02850        if (ifit==nhit-1) tctmp = besttc;
02851        xfit[ifit] = tctmp->GetZPos();
02852        yfit[ifit] = tctmp->GetTPos();
02853      }
02854      LinearFit::Unweighted(nhit,xfit,yfit,parm);
02855      Double_t maxresid = 0.;
02856      for(Int_t ifit=0; ifit<nhit; ++ifit){
02857        maxresid = max(maxresid,
02858                       TMath::Abs(parm[0]+parm[1]*xfit[ifit]-yfit[ifit]));
02859      }
02860      delete [] xfit;
02861      delete [] yfit;
02862      MSG("TrackSR", Msg::kDebug)<< "maxresid = " << maxresid << " / " 
02863                                 << fParmTrk2DMaxResid<< endl;
02864      if (maxresid>fParmTrk2DMaxResid) {
02865        removetc = kFALSE;
02866        planesBefore = 9999999;
02867      }
02868    }//end if not enough hits for anything but a hough track
02869   if(removetc){
02870     // solution found, need to remove bad track cluster(s) between
02871     //the plane of the besttc and the plane represented by
02872     //planesBefore
02873     for(int iremove=0; iremove<planesBefore; ++iremove) {
02874            MSG("TrackSR", Msg::kDebug)<< "removing cluster on plane  " << track->GetCluster(track->GetLast())->GetPlane() << endl;
02875     track->RemoveAt(track->GetLast());
02876       track->Compress();
02877     }
02878   } 
02879   else{
02880     //adding this cluster will somehow make the fit worse, so dont add it
02881     besttc = 0;
02882   }
02883   return removetc;
02884 }

void AlgTrackSRList::CleanUp ( TObjArray *  trackList,
TObjArray *  clusterList 
) [private]

Definition at line 4200 of file AlgTrackSRList.cxx.

References Track2DSR::Compress(), fMapIsWide, TrackClusterSR::GetPlane(), Msg::kDebug, MSG, Track2DSR::RemoveAt(), and tc.

Referenced by RunAlg().

04201 {
04202   Track2DSR *track2d = 0;
04203 
04204   for(Int_t i=0; i<=trackList->GetLast(); ++i){
04205     track2d = dynamic_cast<Track2DSR*>(trackList->At(i));
04206     trackList->RemoveAt(i);
04207     delete track2d;
04208   }
04209   trackList->Compress();
04210   
04211   // delete track clusters
04212   TrackClusterSR *tc = 0;
04213   for(Int_t i=0; i<=clusterList->GetLast(); ++i){
04214    
04215     tc = dynamic_cast<TrackClusterSR*>(clusterList->At(i));
04216     fMapIsWide[tc->GetPlane()] = 0;
04217    
04218     MSG("TrackSR", Msg::kDebug) << "cleaning up, fMapIsWide for plane " 
04219                                << tc->GetPlane() << " is " 
04220                                << fMapIsWide[tc->GetPlane()]
04221                                << endl;
04222     delete tc;
04223   }
04224   clusterList->Compress();
04225   
04226   return;
04227 }

void AlgTrackSRList::DeleteTwinTracks ( CandTrackSRHandleItr &  candTracks1,
CandTrackSRHandleItr &  candTracks2,
TObjArray *  candTracks,
Int_t  uMaxTrack,
Int_t  vMaxTrack,
Int_t  nmax3dtrk,
CandHandle ch 
) [private]

Definition at line 3961 of file AlgTrackSRList.cxx.

References fParmTrk3DTrackMax, fParmTrk3DTwinFrac, CandRecoHandle::GetBegPlane(), Track2DSR::GetChi2(), Track2DSR::GetCluster(), CandRecoHandle::GetEndPlane(), Track2DSR::GetLast(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), CandRecoHandle::GetNPlane(), CandRecoHandle::GetNStrip(), TrackClusterSR::GetNStrip(), TrackClusterSR::GetPlane(), CandTrackSRHandle::GetUTrack(), CandTrackSRHandle::GetVTrack(), Msg::kDebug, len, max, min, MSG, and CandHandle::RemoveDaughter().

Referenced by RunAlg().

03966 {
03967 
03968 //   MSG("TrackSR",Msg::kDebug) << "deleting twin tracks" << endl;;
03969  
03970 
03971   Int_t ncompat=0;
03972   Bool_t * compatlist = new Bool_t [candTracks1.SizeSelect()];
03973 
03974   //initialize the array
03975   Int_t ntrk = candTracks1.SizeSelect();
03976 
03977   for(int i=0;i<ntrk;i++){
03978     compatlist[i]=false;
03979   }
03980   
03981   candTracks1.ResetFirst();
03982   candTracks2.ResetFirst();
03983 
03984   CandTrackSRHandle *track0 = 0;
03985   CandTrackSRHandle *besttrack = 0;
03986   
03987   Int_t itrk0 = 0;
03988 
03989   Int_t besttrk = 0;
03990   Int_t bestlen=0;
03991   Float_t bestchi2=1e9;
03992   Int_t bestnplanes=0;
03993   // find the best track, based on len, and chi2
03994   while( (track0 = candTracks1()) ){
03995     Int_t nplanes=0;
03996     Int_t len=0;
03997     Float_t chi2=0;   
03998     len =TMath::Abs(track0->GetBegPlane()-track0->GetEndPlane());
03999     nplanes = track0->GetNPlane();
04000     chi2 =track0->GetUTrack()->GetChi2()+track0->GetVTrack()->GetChi2();
04001     if(nplanes>bestnplanes){
04002       bestnplanes=nplanes;
04003       bestlen=len;
04004       bestchi2=chi2;
04005       besttrack=track0;
04006       besttrk=itrk0;
04007     }
04008     else if (nplanes==bestnplanes){
04009       if(len>bestlen){
04010         bestnplanes=nplanes;
04011         bestlen=len;
04012         bestchi2=chi2;
04013         besttrack=track0;
04014         besttrk=itrk0;
04015       }
04016       else if( len==bestlen){
04017         if(chi2<bestchi2){
04018           bestnplanes=nplanes;
04019           bestlen=len;
04020           bestchi2=chi2;
04021           besttrack=track0;
04022           besttrk=itrk0;
04023         }
04024       }
04025     }
04026     itrk0++;
04027   }
04028   if(besttrack){
04029     MSG("TrackSR", Msg::kDebug) << "best track  " << besttrk << " extent: " << besttrack->GetEndPlane() << " " << besttrack->GetBegPlane() << endl;
04030   }
04031   if(!besttrack)return;
04032   ncompat=1;
04033   compatlist[besttrk]=true;
04034 
04035   // now find all other tracks compatible with the best
04036   candTracks1.Reset();
04037   itrk0=0;
04038   Int_t itrk1=0;
04039   const CandTrackSRHandle * track1=0;
04040   while( (track0 = candTracks1()) ){
04041     MSG("TrackSR", Msg::kDebug) << "outer loop track " << itrk0 << endl;
04042     Bool_t trackcompat=true;    
04043     for(int itrk=0;itrk<ntrk;itrk++){
04044       itrk1=itrk;
04045       TObject *tobj1 = candTracks->At(itrk);
04046       track1 = dynamic_cast<CandTrackSRHandle*>(tobj1);
04047       if(compatlist[itrk1] && *track0!=*track1){
04048         MSG("TrackSR", Msg::kDebug) << "inner loop track " << itrk1 << endl;
04049         Double_t x0 = 0.;
04050         Double_t x1 = 0.;
04051         Double_t x01 = 0.;
04052         
04053         Double_t x0v = 0.;
04054         Double_t x1v = 0.;
04055         Double_t x01v = 0.;
04056         
04057         Double_t x0u = 0.;
04058         Double_t x1u = 0.;
04059         Double_t x01u = 0.;
04060         x0u=track0->GetUTrack()->GetLast()+1;
04061         x1u=track1->GetUTrack()->GetLast()+1;
04062         TrackClusterSR *clus1;
04063         TrackClusterSR *clus2;
04064         for(Int_t i = 0; i<=track0->GetUTrack()->GetLast(); ++i){
04065           clus1 = track0->GetUTrack()->GetCluster(i);
04066           for(Int_t j = 0; j<=track1->GetUTrack()->GetLast(); ++j){
04067             clus2 = track1->GetUTrack()->GetCluster(j);
04068             if((clus1->GetPlane()==clus2->GetPlane()) && (clus1->GetMinStrip()==clus2->GetMinStrip()) && 
04069                (clus1->GetMaxStrip()==clus2->GetMaxStrip()) && (clus1->GetNStrip()==clus2->GetNStrip())) {
04070               x01u++;
04071             }
04072           }      
04073         }        
04074         
04075         x0v=track0->GetVTrack()->GetLast()+1;
04076         x1v=track1->GetVTrack()->GetLast()+1;
04077         for(Int_t i = 0; i<=track0->GetVTrack()->GetLast(); ++i){
04078           clus1 = track0->GetVTrack()->GetCluster(i);
04079           for(Int_t j = 0; j<=track1->GetVTrack()->GetLast(); ++j){
04080             clus2 = track1->GetVTrack()->GetCluster(j);
04081             if((clus1->GetPlane()==clus2->GetPlane()) && (clus1->GetMinStrip()==clus2->GetMinStrip()) && 
04082                (clus1->GetMaxStrip()==clus2->GetMaxStrip()) && (clus1->GetNStrip()==clus2->GetNStrip())) {              
04083               x01v++;
04084             }
04085           }      
04086         }   
04087         x0 = 1.*track0->GetNStrip();
04088         x1 = 1.*track1->GetNStrip();
04089         x01 = 1.*track0->GetNStrip(track1);                       
04090         
04091         MSG("TrackSR", Msg::kDebug) << " V overlap " << x01v/x0v << "/" << fParmTrk3DTwinFrac << "/" << x01v/x1v << " U overlap " << x01u/x0u << "/" << fParmTrk3DTwinFrac << "/" << x01u/x1u << endl;  
04092         MSG("TrackSR", Msg::kDebug) << " 3d  overlap " << x01/x0 << "/" << fParmTrk3DTwinFrac << "/" << x01/x1  << endl;
04093         if(x01v/x0v>=fParmTrk3DTwinFrac ||
04094            x01v/x1v>=fParmTrk3DTwinFrac ||
04095            x01u/x0u>=fParmTrk3DTwinFrac ||
04096            x01u/x1u>=fParmTrk3DTwinFrac ||
04097            x01/x0>=fParmTrk3DTwinFrac  ||
04098            x01/x1>=fParmTrk3DTwinFrac) {
04099           MSG("TrackSR", Msg::kDebug) << " tracks have too much overlap - eliminate track  " << itrk0 << endl;
04100           trackcompat=false;
04101           Int_t len0 =TMath::Abs(track0->GetBegPlane()-track0->GetEndPlane());
04102           Int_t nplanes0 = track0->GetNPlane();
04103           Float_t chi2_0 =track0->GetUTrack()->GetChi2()+track0->GetVTrack()->GetChi2();
04104           
04105           Int_t len1 = TMath::Abs(track1->GetBegPlane()-track1->GetEndPlane());
04106           Int_t nplanes1 = track1->GetNPlane();
04107           Float_t chi2_1 =track1->GetUTrack()->GetChi2()+track1->GetVTrack()->GetChi2();
04108           
04109           if(nplanes0>nplanes1){
04110             trackcompat=true;
04111             compatlist[itrk0]=true;
04112             compatlist[itrk1]=false;
04113             MSG("TrackSR", Msg::kDebug) << " eliminating track  " << itrk1 << " instead " << endl;
04114           }
04115           else if (nplanes0==nplanes1){
04116             if(len0>len1){
04117               trackcompat=true;
04118               compatlist[itrk0]=true;
04119               compatlist[itrk1]=false;
04120               MSG("TrackSR", Msg::kDebug) << " eliminating track  " << itrk1 << "instead " << endl;
04121             }
04122             else if( len0==len1){
04123               if(chi2_0<chi2_1){
04124                 trackcompat=true;
04125                 compatlist[itrk0]=true;
04126                 compatlist[itrk1]=false;
04127                 MSG("TrackSR", Msg::kDebug) << " eliminating track " << itrk1 << "instead " << endl;
04128               }
04129             }
04130           }
04131           if(!compatlist[itrk0])break;
04132         }
04133       }
04134     }
04135     if(trackcompat && !compatlist[itrk0])ncompat++;
04136     compatlist[itrk0]=trackcompat;
04137     itrk0++;
04138   }
04139   MSG("TrackSR", Msg::kDebug) << ncompat << " compatible tracks " << endl;
04140   // sort tracks by strip multiplicity
04141   Int_t *iassigned = new Int_t[candTracks->GetLast()+1];
04142   for (Int_t indx=0; indx<=candTracks->GetLast(); ++indx) {
04143     iassigned[indx] = 0;
04144   }
04145 
04146   Int_t nmax3d = 0;
04147   if(fParmTrk3DTrackMax>0) nmax3d = fParmTrk3DTrackMax;
04148   else{
04149     nmax3d = max(uMaxTrack,vMaxTrack);
04150     nmax3d = max(nmax3d,nmax3dtrk);
04151   }
04152 
04153   // allow for at least one track
04154   if (nmax3d<1) nmax3d = 1;
04155 
04156   MSG("TrackSR", Msg::kDebug) << "nmax3d = " << nmax3d <<  " p3DTrackMax  " << fParmTrk3DTrackMax << endl;
04157 
04158   //loop over all the tracks to find the n tracks with the most
04159   //hit strips where n = minimum(nmax3d, candTracks->GetLast()+1
04160   Int_t nLoop = min(nmax3d,candTracks->GetLast()+1);
04161   for(Int_t i = 0; i < nLoop; ++i){
04162     int imax=-1;
04163     Int_t nmax = 0;
04164     
04165     //loop over the tracks that werent removed to find which one has 
04166     //the largest number of strips - that will be the one we eventually keep
04167     for (itrk0=0; itrk0<=candTracks->GetLast(); ++itrk0){
04168       TObject *tobj0 = candTracks->At(itrk0);
04169       const CandTrackSRHandle *track0 =
04170         dynamic_cast<CandTrackSRHandle*>(tobj0);
04171       if (compatlist[itrk0] && !iassigned[itrk0]){
04172         if (track0->GetNStrip()>nmax) {
04173           nmax = track0->GetNStrip();
04174           imax = itrk0;
04175         }
04176       }
04177     }//end loop over tracks
04178     if(imax>=0) iassigned[imax] = 1;
04179   }//end loop to get n tracks with most hit strips
04180 
04181   //get rid of all tracks but the one with the most strips hit in it
04182   for(itrk0=0; itrk0<=candTracks->GetLast(); ++itrk0){
04183     if(!iassigned[itrk0]){
04184       // remove from candidate daughter list
04185       MSG("TrackSR", Msg::kDebug) << "removing track at position " << itrk0+1 
04186                                   << " of "
04187                                  << candTracks->GetLast()+1 << endl;
04188       ch.RemoveDaughter(dynamic_cast<CandTrackSRHandle*>(candTracks->At(itrk0)));
04189     }
04190   }
04191   
04192   delete [] iassigned;
04193   delete [] compatlist;
04194   return;
04195 }

void AlgTrackSRList::FillInGaps ( Track2DSRItr &  trackItr,
TrackClusterSRItr &  clusterItr,
Int_t  direction 
) [private]

Definition at line 3525 of file AlgTrackSRList.cxx.

References Track2DSR::Add(), fParmTrk2DMaxResid, fParmTrk2DNSkip, Track2DSR::GetBegPlane(), TrackClusterSR::GetCharge(), Track2DSR::GetCluster(), Track2DSR::GetEndPlane(), Track2DSR::GetLast(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), Track2DSR::GetPlaneView(), TrackClusterSR::GetTPos(), TrackClusterSR::GetTPosError(), TrackClusterSR::GetZPos(), Msg::kDebug, PlaneView::kUnknown, min, MSG, Track2DSR::RemoveAt(), and tc.

Referenced by RunAlg().

03528 {
03529   Track2DSR *thistrack = 0;
03530   MSG("TrackSR",Msg::kDebug) << "filling in gaps" << endl;
03531 
03532   //array of clusters to add to the track
03533   TObjArray *addtotrack = new TObjArray(1,0);
03534 
03535   TrackClusterSR *tc = 0;
03536   TrackClusterSR *tc0 = 0;
03537   TrackClusterSR *tc1 = 0;
03538   TrackClusterSR *besttc = 0;
03539 
03540   clusterItr.ResetFirst();
03541   
03542   PlaneView::PlaneView_t view = PlaneView::kUnknown;
03543   if(clusterItr.Ptr() ) view = clusterItr.Ptr()->GetPlaneView();
03544 
03545   Int_t numAdded = 0;
03546 
03547   Double_t dslope = 0.;
03548   Double_t tcerr0 = 0.;
03549   Double_t tcerr1 = 0.;
03550   Double_t tcerr = 0.;
03551   Double_t bestresid = 99999999.;
03552   Double_t proj = 0.;
03553   Double_t dtpos = 0.;
03554  
03555   trackItr.ResetFirst();
03556   clusterItr.Reset();
03557   
03558 
03559   if(direction == 1)clusterItr.ResetFirst();
03560   else if(direction == -1) clusterItr.ResetLast(); 
03561   trackItr.ResetFirst();
03562 
03563   //loop over tracks
03564   while( (thistrack = trackItr()) ){
03565 
03566     //on a new track so clear the addtotrack array
03567     addtotrack->Clear();
03568     numAdded = 1;
03569 
03570     map<TrackClusterSR*,Double_t> thisslope;
03571 
03572     while( numAdded>0 && thistrack->GetPlaneView() == view){
03573       
03574       //reset numAdded to 0 for the loop
03575       numAdded = 0;
03576       
03577  
03578         
03579       MSG("TrackSR",Msg::kDebug) << "track beg plane = " 
03580                                         << thistrack->GetBegPlane() 
03581                                         << " end plane = " 
03582                                         << thistrack->GetEndPlane() 
03583                                  << endl;
03584       
03585       //loop over the clusters in this track
03586       for(int ihit=0; ihit<thistrack->GetLast(); ++ihit){
03587         tc0 = thistrack->GetCluster(ihit);
03588         tc1 = thistrack->GetCluster(ihit+1);
03589         
03590         //make sure clusters are not from consecutive planes in the same view
03591         //also if one of the clusters is a track end point, look for clusters
03592         //beyond it
03593         if(TMath::Abs(tc0->GetPlane()-tc1->GetPlane())>2
03594            || ihit == 0 || ihit+1 == thistrack->GetLast() ){
03595           
03596           MSG("TrackSR",Msg::kDebug) << "tc0 " << tc0->GetPlane() << " " 
03597                                             << tc0->GetMinStrip() << " " 
03598                                             << tc0->GetMaxStrip() << " " 
03599                                             << tc0->GetTPos() << " " 
03600                                             << tc0->GetTPosError() << endl 
03601                                             << "tc1 " << tc1->GetPlane() << " " 
03602                                             << tc1->GetMinStrip() << " " 
03603                                             << tc1->GetMaxStrip() << " " 
03604                                             << tc1->GetTPos() << " " 
03605                                             << tc1->GetTPosError() << endl;
03606           
03607           dslope = (tc0->GetTPos()-tc1->GetTPos())/(tc0->GetZPos()-tc1->GetZPos());
03608           tcerr0 = tc0->GetTPosError();
03609           tcerr1 = tc1->GetTPosError();
03610           tcerr = sqrt(tcerr0*tcerr0+tcerr1*tcerr1);
03611           
03612           // limit tcerr to 5 cm
03613           if (tcerr>0.05) tcerr = 0.05;
03614           
03615           MSG("TrackSR",Msg::kDebug) << "dslope " << dslope << "  errors = " 
03616                                      << tcerr0 << " " << tcerr1 << " " << tcerr 
03617                                      << endl;
03618           
03619           //on a new gap, so reset the best residual and cluster variables
03620           bestresid = 99999.;
03621           besttc = 0;
03622           
03623           //loop over the master list of track clusters
03624           if(direction == 1)clusterItr.ResetFirst();
03625           else if(direction == -1) clusterItr.ResetLast();
03626           
03627           while( (tc = clusterItr()) ){
03628 
03629             
03630             //check to see if tc comes from the plane before the first in 
03631             //the track and isnt too far away from the first hit plane
03632             
03633             //or if tc comes from the plane after the last in
03634             //the track and isnt too far away from the last hit plane
03635             
03636             //or if tc comes from a plane between tc0 and tc1
03637             
03638             //remember, the track's clusters are in order according to the
03639             //direction, so the plane of tc0 is before that of tc1 in time
03640 
03641 
03642             if((tc->GetPlane()*direction>tc0->GetPlane()*direction
03643                 && tc->GetPlane()*direction<tc1->GetPlane()*direction)
03644                  || (ihit+1 == thistrack->GetLast() 
03645                      && tc->GetPlane()*direction>tc1->GetPlane()*direction
03646                      && TMath::Abs(tc->GetPlane()
03647                                    -tc1->GetPlane())<=fParmTrk2DNSkip*TMath::Abs(tc0->GetPlane()-tc1->GetPlane()))
03648                  || (ihit == 0 && 
03649                      tc->GetPlane()*direction<tc0->GetPlane()*direction
03650                      && TMath::Abs(tc->GetPlane()
03651                                    -tc0->GetPlane())<=fParmTrk2DNSkip*TMath::Abs(tc0->GetPlane()-tc1->GetPlane()))){
03652               
03653                 
03654                 //get the projected location for the plane of tc from the slope
03655                 //between the two bounding clusters
03656                 
03657                 proj = tc0->GetTPos()+dslope*(tc->GetZPos()-tc0->GetZPos());
03658                 
03659                 dtpos = 0.;
03660                 
03661                 //if the transverse position of tc is out of bounds based on the
03662                 //projected location +- the error, then set dtpos to the minimum
03663                 //of the difference in its tpos +- its error - projection +- error
03664                 if(tc->GetTPos()<proj-tcerr || tc->GetTPos()>proj+tcerr)
03665                   dtpos=min(TMath::Abs(tc->GetTPos()+tc->GetTPosError()
03666                                        -(proj+tcerr)),
03667                             TMath::Abs(tc->GetTPos()-tc->GetTPosError()
03668                                        -(proj-tcerr)));
03669                 
03670                 MSG("TrackSR",Msg::kDebug) << "considering " 
03671                                                   << tc->GetPlane() 
03672                                                   << " tpos " << tc->GetTPos() 
03673                                                   << " proj " << proj 
03674                                                   << " dtpos " << dtpos 
03675                                            << "/" <<  fParmTrk2DMaxResid*(1+0.15*TMath::Abs(tc->GetPlane()-tc0->GetPlane()))
03676                                                   << endl;
03677                 
03678                 //if the dtpos is acceptable and less than the best residual so far
03679                 Float_t dplaneproj=min(TMath::Abs(tc->GetPlane()-tc0->GetPlane()),TMath::Abs(tc->GetPlane()-tc1->GetPlane()));
03680                 if(dtpos<fParmTrk2DMaxResid*(1+0.15*dplaneproj) && dtpos<bestresid){
03681                   besttc = tc;
03682                   bestresid = dtpos;
03683                 }
03684                 //            }
03685             }//end if on a plane between the bounding clusters
03686           }//end loop over master cluster list
03687           
03688           //you found a cluster to add, so add it
03689           if(besttc){
03690             addtotrack->AddLast(besttc);
03691             thisslope[besttc] = ((tc0->GetTPos()-tc1->GetTPos())
03692                                  /(tc0->GetZPos()-tc1->GetZPos()));
03693             ++numAdded;
03694             MSG("TrackSR",Msg::kDebug) << "adding " << besttc->GetPlane() 
03695                                               << " " << besttc->GetTPos() 
03696                                               << " numAdded = " << numAdded <<endl;
03697           }
03698         }//end if in a gap between clusters
03699       }//end loop over clusters in the track
03700 
03701       
03702       //now that you have found all the clusters to add to the track, add them
03703       for(int i=0; i<=addtotrack->GetLast(); ++i) {
03704         tc1 = dynamic_cast<TrackClusterSR*>(addtotrack->At(i));
03705         MSG("TrackSR",Msg::kDebug) << " adding " << tc1->GetPlane() 
03706                                           << " " << tc1->GetTPos() << endl;
03707         thistrack->Add(tc1,thisslope[tc1]);
03708 
03709         //added the cluster, now remove it from the addtotrack array
03710         addtotrack->RemoveAt(i);
03711       }
03712 
03713       addtotrack->Compress();
03714     
03715       }//end loop to add clusters to the track
03716 
03717 
03718       // remove isloated clusters at track ends i
03719     if(thistrack->GetLast()>2){
03720       Float_t dz_front=TMath::Abs(thistrack->GetCluster(0)->GetZPos()-
03721                             thistrack->GetCluster(1)->GetZPos());
03722       while(dz_front>3.0 && thistrack->GetLast()>2){
03723         thistrack->RemoveAt(0);
03724         dz_front=TMath::Abs(thistrack->GetCluster(0)->GetZPos()-
03725                       thistrack->GetCluster(1)->GetZPos());
03726       }
03727     }
03728     if(thistrack->GetLast()>2){
03729       Float_t dz_back=TMath::Abs(thistrack->GetCluster(thistrack->GetLast())->GetZPos()-
03730                            thistrack->GetCluster(thistrack->GetLast()-1)->GetZPos());
03731       while(dz_back>3.0 && thistrack->GetLast()>2){
03732         thistrack->RemoveAt(thistrack->GetLast());
03733         dz_back=TMath::Abs(thistrack->GetCluster(thistrack->GetLast())->GetZPos()-
03734                      thistrack->GetCluster(thistrack->GetLast()-1)->GetZPos());    
03735       }
03736     }
03737     Float_t totcharge = 0.;
03738     for (int i=0; i<=thistrack->GetLast(); i++) {
03739       tc1 = thistrack->GetCluster(i);
03740       MSG("TrackSR",Msg::kDebug) << " " << tc1->GetPlane() 
03741                                         << " " << tc1->GetTPos() << endl;
03742       totcharge += tc1->GetCharge();
03743     }
03744     //     MSG("TrackSR",Msg::kVerbose) << "total charge = " << totcharge << "\n";
03745     
03746     }//end loop over tracks
03747 
03748   delete addtotrack;
03749   return;
03750 }

Int_t AlgTrackSRList::FillTrackList ( TrackClusterSRItr &  clusterItr,
TObjArray *  trackList,
TObjArray *  seedClusters,
Double_t  houghSlope,
Double_t  houghIntercept,
Int_t  direction,
Int_t  iterate,
Int_t &  begPlane 
) [private]

Definition at line 1960 of file AlgTrackSRList.cxx.

References Track2DSR::Add(), Track2DSR::Compress(), Track2DSR::fIterate, fParmHitNTime, fParmHitTime, fParmSingleHitDef, fParmTrk2DPlnEnd, fParmTrk2DWin0, TrackClusterSR::GetBegTime(), TrackClusterSR::GetCharge(), TrackClusterSR::GetNStrip(), TrackClusterSR::GetPlane(), TrackClusterSR::GetTPos(), TrackClusterSR::GetZPos(), TrackClusterSR::IsValid(), TrackClusterSR::IsWide(), Msg::kDebug, min, MSG, Track2DSR::RemoveAt(), Track2DSR::SetDirection(), Track2DSR::SetHoughIntercept(), Track2DSR::SetHoughSlope(), TrackClusterSR::SetIsValid(), and tc.

Referenced by RunAlg().

01967 {
01968 
01969   TrackClusterSR *cluster = 0;
01970   TrackClusterSR *nextCluster=0;
01971   TrackClusterSRItr nextclusterItr2(clusterItr);
01972 
01973   Int_t prevTrackCount = trackList->GetLast()+1;
01974   MSG("TrackSR", Msg::kDebug) << "previously " << prevTrackCount 
01975                              << " tracks in list" << endl;
01976  //set the direction of the iterator based on the timing direction.
01977   //then the first plane represent by clusters is your first plane
01978   //that is a possible seed hit.
01979   if(direction == 1) clusterItr.ResetFirst();
01980   else if(direction == -1) clusterItr.ResetLast();
01981 
01982   vector<Int_t> planeindex;
01983   Int_t ncount=0;
01984   Int_t oldipln=0;
01985   while (TrackClusterSR *tc = clusterItr()) {
01986     Int_t ipln = tc->GetPlane();
01987 
01988     if(tc->IsValid()){
01989       if (!ncount || ipln!=oldipln ) {
01990         planeindex.push_back(ipln);
01991         MSG("TrackSR", Msg::kDebug) << " valid cluster on plane " << ipln << endl;
01992       }
01993       ncount++;
01994       oldipln = ipln;
01995     }
01996   }
01997   
01998   //difference in plane number between plane of current cluster
01999   //and one being looked at
02000   Int_t deltaPlane = 0;
02001 
02002   //difference in time between current cluster and one being looked at
02003   Double_t deltaTime = 0.;
02004   
02005   //difference in z pos and transverse position between current cluster
02006   //and one being looked at
02007   Double_t deltaTPos = 0.;
02008   Double_t deltaZPos = 0.;
02009 
02010   //keep track of which valid cluster you are on
02011   Int_t validCtr = 0;
02012   //loop over all the clusters in the list.  only those that are still
02013   //valid will be considered to be used as new seed clusters.  non-valid
02014   //clusters have alread been used in a track
02015   for (Int_t iplnindx=0; iplnindx<static_cast<Int_t>(planeindex.size());iplnindx++) {
02016     Int_t ipln = planeindex[iplnindx];
02017     MSG("TrackSR",Msg::kDebug) << "PLANE " << ipln  << "\n";
02018     //reset the direction of the iterator.
02019     if(direction == 1) clusterItr.ResetFirst();
02020     else if(direction == -1) clusterItr.ResetLast();
02021     while( (cluster = clusterItr()) ){
02022       if(cluster->GetPlane()==ipln && cluster->IsValid()){
02023         //if you have a valid cluster, it could be a seed hit
02024          
02025         MSG("TrackSR",Msg::kDebug) << cluster->GetPlane() 
02026                                    << " " << cluster->GetTPos() << " " 
02027                                    << cluster->GetZPos()
02028                                    << " " << cluster->GetCharge() << " " 
02029                                    << (Int_t)cluster->IsValid() << " " 
02030                                    << cluster->GetNStrip() << endl;
02031         
02032         //seed clusters are those which are well separated from clusters
02033         //preceding them in the list in terms of plane number.  they must 
02034         //also be well separated in time and transverse position
02035         
02036         if(iplnindx==0){
02037           ++validCtr;
02038           //first valid plane in the bunch, call it begPlane
02039           begPlane = cluster->GetPlane();
02040           //first valid cluster in the view, it is a seed cluster
02041           Track2DSR *track = new Track2DSR();
02042           MSG("TrackSR",Msg::kDebug) << " Found Seed hit " << cluster->GetNStrip() << "/" << fParmSingleHitDef ;
02043           //set the fIterate variable for track so it can id
02044           //the iteration it was created in
02045           track->fIterate = iterate;
02046           track->SetDirection(direction);
02047           if(!cluster->IsWide()){
02048             MSG("TrackSR",Msg::kDebug) << " OK as seed" << endl;
02049             seedClusters->AddLast(cluster);
02050           }
02051           //too many hits in the cluster to use as a seed cluster
02052           //set the cluster as non-valid so it doesnt get used
02053           //in the initial finding somehow.  later when we loop over
02054           //all tracks to fill in the gaps it could be used
02055           else{
02056             MSG("TrackSR",Msg::kDebug) << " too wide for seed " << endl;
02057             track->fIterate = -1;
02058             cluster->SetIsValid(false);
02059           }
02060           track->Add(cluster,houghSlope);
02061           track->SetHoughSlope(houghSlope);
02062           track->SetHoughIntercept(houghIntercept);
02063           trackList->AddLast(track);
02064         }
02065         else{
02066           Bool_t isolated=true;
02067     //reset the direction of the iterator.
02068           if(direction == 1) nextclusterItr2.ResetFirst();
02069           else if(direction == -1) nextclusterItr2.ResetLast();
02070           while( (nextCluster = nextclusterItr2()) ){
02071             if(nextCluster->GetPlane()==begPlane && nextCluster->IsValid()){
02072               deltaPlane = TMath::Abs(cluster->GetPlane() - nextCluster->GetPlane());
02073               deltaTime = TMath::Abs(cluster->GetBegTime() - nextCluster->GetBegTime());
02074               deltaZPos = TMath::Abs(cluster->GetZPos() - nextCluster->GetZPos());
02075               deltaTPos = min(TMath::Abs(cluster->GetTPos() - nextCluster->GetTPos() 
02076                                      + deltaZPos*houghSlope),TMath::Abs(cluster->GetTPos() - nextCluster->GetTPos()));
02077               
02078               MSG("TrackSR", Msg::kDebug) << "current cluster plane = " 
02079                 
02080                                           << cluster->GetPlane()
02081                                           << " prev cluster plane = " 
02082                                           << nextCluster->GetPlane()
02083                                           << endl << " deltaPlane = " << deltaPlane 
02084                                           << " / " << fParmTrk2DPlnEnd << endl  
02085                                           << " deltaTime = " << deltaTime << "/ " 
02086                                           << fParmHitNTime << "/" << fParmHitTime
02087                                           << " deltaTPos = " << deltaTPos << " / " 
02088                                           << fParmTrk2DWin0*deltaZPos << endl;
02089               
02090               if(deltaPlane<=fParmTrk2DPlnEnd 
02091                  && deltaTime>fParmHitNTime && deltaTime<fParmHitTime
02092                  && deltaTPos<fParmTrk2DWin0*deltaZPos)
02093                 {
02094                   isolated=false;
02095                   break;
02096                 }
02097             }
02098           }
02099           if(isolated){
02100             //make a new track
02101             Track2DSR *track = new Track2DSR();
02102             MSG("TrackSR",Msg::kDebug) << " Found Seed hit " << endl;
02103             //set the fIterate variable for track so it can id
02104             //the iteration it was created in
02105             track->fIterate = iterate;
02106             track->SetDirection(direction);
02107             if(!cluster->IsWide()){
02108               seedClusters->AddLast(cluster);
02109             }
02110             //too many hits in the cluster to use as a seed cluster
02111             //set the cluster as non-valid so it doesnt get used
02112             //in the initial finding somehow.  later when we loop over
02113             //all tracks to fill in the gaps it could be used
02114             else{
02115               track->fIterate = -1;
02116               cluster->SetIsValid(false);
02117             }
02118             track->Add(cluster,houghSlope);
02119             track->SetHoughSlope(houghSlope);
02120             track->SetHoughIntercept(houghIntercept);
02121             trackList->AddLast(track);
02122           }       
02123         }
02124       }
02125     }
02126     begPlane=ipln;
02127   }
02128 
02129   //loop over the tracks one more time to get rid of the flagged tracks
02130   //i could have maybe did this in the above loop, but its not clear
02131   //that i wouldnt end up deleting an object before the last time it would
02132   //be accessed by the loop
02133 
02134   Track2DSR *track1 = 0;
02135   for(Int_t i = 0; i<=trackList->GetLast(); ++i){
02136     track1 = dynamic_cast<Track2DSR *>(trackList->At(i));
02137     if(track1->fIterate < 0){
02138       trackList->RemoveAt(i);
02139       delete track1;
02140     }
02141   }//end loop to remove tracks
02142   
02143   trackList->Compress();
02144   return trackList->GetLast()-prevTrackCount;
02145 }

void AlgTrackSRList::Find2DTrackEndPoints ( CandStripHandleItr &  stripItr,
Int_t &  vtxPlane,
Int_t &  endPlane,
Int_t &  numPlanes 
) [private]

Definition at line 4233 of file AlgTrackSRList.cxx.

References fParmMinStripPulseHeight, CandStripHandle::GetCharge(), CandStripHandle::GetPlane(), CandStripHandle::GetTime(), CandStripHandle::GetTPos(), CandStripHandle::GetZPos(), Msg::kDebug, and MSG.

Referenced by RunAlg().

04236 {
04237   //prevPlane keeps track of which plane you were last on
04238   Int_t prevPlane = 0;
04239   CandStripHandle *strip = 0;
04240   Int_t currentPlane = 0;
04241 
04242   while( stripItr.IsValid() ){
04243     
04244     strip = stripItr.Ptr();
04245     currentPlane = strip->GetPlane();
04246 
04247     MSG("TrackSR", Msg::kDebug) << "plane " << currentPlane << " tpos " 
04248                                << strip->GetTPos() << " zpos " << strip->GetZPos()
04249                                << " time " << strip->GetTime() << endl;
04250 
04251     if( prevPlane != currentPlane && 
04252         strip->GetCharge()>fParmMinStripPulseHeight){
04253       ++numPlanes;
04254       if( currentPlane < vtxPlane) vtxPlane = currentPlane;
04255       if( currentPlane > endPlane) endPlane = currentPlane;
04256       prevPlane = currentPlane;
04257     }//end if this is a new plane
04258     
04259     //advance the iterator
04260     stripItr.Next();
04261   }//end loop over u strips
04262   
04263   stripItr.ResetFirst();
04264 
04265   return;
04266 }

Int_t AlgTrackSRList::FindNumSkippedPlanes ( Int_t  currentPlane,
Int_t  prevPlane,
Float_t  PrevUPos,
Float_t  PrevVPos,
Float_t  prevZPos,
Double_t  trackSlopeU,
Double_t  trackSlopeV,
Int_t  direction,
PlaneView::PlaneView_t  view 
) [private]

Definition at line 2522 of file AlgTrackSRList.cxx.

References fDetector, fMapIsWide, fParmExtrapError, fParmTrk2DNSkip, fVldc, PlexPlaneId::GetAdjoinScint(), PlexPlaneId::GetPlane(), PlexPlaneId::GetPlaneCoverage(), PlexPlaneId::GetPlaneView(), UgliGeomHandle::GetScintPlnHandle(), UgliPlnHandle::GetZ0(), PlexPlaneId::IsValid(), UgliScintPlnHandle::IsValid(), Msg::kDebug, MSG, and PlaneIsActive().

02529 { 
02530   UgliGeomHandle ugh(fVldc);
02531 
02532   Int_t thisPlane = prevPlane+direction;
02533   Int_t numSkipped = 0;
02534 
02535   //projected transverse position at the currentPlane
02536   //  Float_t tPosProj = 999.;
02537 
02538   Float_t uPosProj = prevUPos;
02539   Float_t vPosProj = prevVPos;
02540   UgliScintPlnHandle planeid;
02541   PlexPlaneId plexPlane(fDetector,prevPlane, 0);
02542   PlexPlaneId plexPlaneEnd(fDetector,currentPlane, 0);
02543   plexPlane = plexPlane.GetAdjoinScint(direction);
02544 
02545   // break out of loop when count exceeds 5 times value of nskip
02546   // parameter (max # skipped planes allowed in spectrometer)
02547 
02548   while(plexPlane != plexPlaneEnd && numSkipped<=fParmTrk2DNSkip*5){
02549     
02550     thisPlane = plexPlane.GetPlane();
02551     //increment numSkipped if any of the above reasons are true
02552     if(plexPlane.IsValid() &&  fMapIsWide[thisPlane]==0 
02553        && plexPlane.GetPlaneView()==view) {
02554       MSG("TrackSR",Msg::kDebug) << "looking at plane " << thisPlane 
02555                                  << " view = " << plexPlane.GetPlaneView() 
02556                                  << " coverage " << plexPlane.GetPlaneCoverage() << endl;
02557       
02558       // !! Consider as "spectrometer planes" the ones that are 
02559       // in the partially instrumented region of the detector, in order 
02560       // to better accomodate track finding there
02561       
02562         planeid = ugh.GetScintPlnHandle(plexPlane);
02563         uPosProj = prevUPos;
02564         vPosProj = prevVPos;
02565         if(planeid.IsValid()){ 
02566           uPosProj = prevUPos + (planeid.GetZ0()-prevZPos)*trackSlopeU;
02567           vPosProj = prevVPos + (planeid.GetZ0()-prevZPos)*trackSlopeV;
02568           
02569             MSG("TrackSR",Msg::kDebug) << " projected U position " << uPosProj 
02570                              << " projected V position " << vPosProj 
02571                              << " prev U position " << prevUPos
02572                              << " prev V position " << prevVPos
02573                              << " iswide plane = " 
02574                              << fMapIsWide[thisPlane] 
02575                              << "  dplaneoff = " << numSkipped << endl 
02576                              << "  slope = " << trackSlopeU << "/" 
02577                              << trackSlopeV 
02578                              << endl;
02579           
02580           Float_t posErr = TMath::Abs(planeid.GetZ0()-prevZPos)*fParmExtrapError;
02581           if(posErr>0.2)posErr=0.2;
02582           if( PlaneIsActive(thisPlane,uPosProj,vPosProj, posErr)){
02583             ++numSkipped;
02584             MSG("TrackSR",Msg::kDebug) << "plane is active " << numSkipped << endl;
02585           }
02586         } // end if planeid is valid 
02587       }//end if plexplane is valid
02588     
02589     plexPlane = plexPlane.GetAdjoinScint(direction);
02590   }//end while looking for missing planes
02591   
02592   return numSkipped;
02593 }

Int_t AlgTrackSRList::FindNumSkippedPlanes ( Int_t  currentPlane,
Int_t  prevPlane,
Float_t  prevTPos,
Float_t  prevZPos,
Double_t  trackSlope,
Int_t  direction,
PlaneView::PlaneView_t  view 
) [private]

Definition at line 2441 of file AlgTrackSRList.cxx.

References fDetector, fMapIsWide, fParmTrk2DNSkip, fVldc, PlexPlaneId::GetAdjoinScint(), PlexPlaneId::GetPlane(), PlexPlaneId::GetPlaneCoverage(), PlexPlaneId::GetPlaneView(), UgliGeomHandle::GetScintPlnHandle(), UgliPlnHandle::GetZ0(), PlexPlaneId::IsValid(), UgliScintPlnHandle::IsValid(), Msg::kDebug, Detector::kFar, Detector::kNear, PlaneCoverage::kNearFull, PlaneCoverage::kNearPartial, PlaneCoverage::kNoActive, PlaneView::kU, PlaneView::kV, and MSG.

02445 { 
02446   UgliGeomHandle ugh(fVldc);
02447 
02448   PlaneView::PlaneView_t view0 = PlaneView::kU;
02449   PlaneView::PlaneView_t view1 = PlaneView::kV;
02450 
02451   Int_t thisPlane = prevPlane+direction;
02452   Int_t numSkipped = 0;
02453 
02454   //projected transverse position at the currentPlane
02455   //  Float_t tPosProj = 999.;
02456   Float_t tPosProj = prevTPos;
02457   UgliScintPlnHandle planeid;
02458   PlexPlaneId plexPlane(fDetector,prevPlane, 0);
02459   PlexPlaneId plexPlaneEnd(fDetector,currentPlane, 0);
02460 
02461   plexPlane = plexPlane.GetAdjoinScint(direction);
02462 
02463   while(plexPlane != plexPlaneEnd && numSkipped<=fParmTrk2DNSkip){
02464     
02465     thisPlane = plexPlane.GetPlane();
02466 
02467     MSG("TrackSR",Msg::kDebug) << "looking at plane " << thisPlane 
02468                               << " view = " << plexPlane.GetPlaneView() 
02469                               << endl;
02470 
02471     // !! Consider as "spectrometer planes" the ones that are 
02472     // in the partially instrumented region of the detector, in order 
02473     // to better accomodate track finding there
02474 
02475     if(plexPlane.IsValid()  
02476        && ( (plexPlane.GetPlaneCoverage() != PlaneCoverage::kNoActive 
02477              &&  fDetector==Detector::kFar) 
02478             || (plexPlane.GetPlaneCoverage()==PlaneCoverage::kNearPartial 
02479                 && ((plexPlane.GetPlaneView()==view0 && (tPosProj > 0. && tPosProj<2.4)) 
02480                     || (plexPlane.GetPlaneView()==view1 && (tPosProj  < 0. && tPosProj>-2.4)))) 
02481             || (plexPlane.GetPlaneCoverage()==PlaneCoverage::kNearFull)) 
02482        && plexPlane.GetPlaneView() == view)     
02483       {           
02484         
02485         planeid = ugh.GetScintPlnHandle(plexPlane);
02486         
02487         tPosProj = prevTPos;
02488         //      tPosProj = 999.;
02489         
02490         if(planeid.IsValid()) 
02491           tPosProj = prevTPos + (planeid.GetZ0()-prevZPos)*trackSlope;
02492         
02493         MSG("TrackSR",Msg::kDebug) << "projected position " << tPosProj 
02494                                    << " iswide plane = " 
02495                                    << fMapIsWide[thisPlane] 
02496                                    << "  dplaneoff = " << numSkipped 
02497                                    << " tpos = " << prevTPos 
02498                                    << "  slope = " << trackSlope 
02499                                    << endl;
02500       
02501         //increment numSkipped if any of the above reasons are true
02502         if(planeid.IsValid() &&  fMapIsWide[thisPlane]==0 
02503            && ((fDetector==Detector::kFar
02504                 && TMath::Abs(tPosProj)>0.3) 
02505                || (fDetector==Detector::kNear 
02506                    && TMath::Abs(tPosProj)>0.8))
02507            ) 
02508           ++numSkipped;
02509         
02510       }//end if plexplane is valid
02511     
02512     plexPlane = plexPlane.GetAdjoinScint(direction);
02513   }//end while looking for missing planes
02514 
02515   return numSkipped;
02516 }

Int_t AlgTrackSRList::FindNumSkippedPlanes ( Int_t  currentPlane,
Int_t  prevPlane,
std::vector< Int_t > &  hitPlanes,
Int_t  direction 
) [private]

Definition at line 2597 of file AlgTrackSRList.cxx.

Referenced by AddClustersToTracks(), MergeTracks(), and SpectrometerTracking().

02600 { 
02601   Int_t numSkipped = 0;
02602 
02603   for(Int_t i = 0; i < static_cast<Int_t>(hitPlanes.size()); ++i){
02604     
02605     if(direction*hitPlanes[i]>direction*prevPlane
02606        && direction*hitPlanes[i]<direction*currentPlane) ++numSkipped;
02607   }
02608 
02609   return numSkipped;
02610 }

Int_t AlgTrackSRList::FindTimingDirection ( TrackClusterSRItr &  clusterItr  )  [private]

Definition at line 1857 of file AlgTrackSRList.cxx.

References fParmMaxTimingResid, TrackClusterSR::GetBegTime(), TrackClusterSR::GetCharge(), TrackClusterSR::GetZPos(), Msg::kDebug, min, MSG, Munits::ns, ArrivalTime::Weight(), and LinearFit::Weighted().

01858 {
01859   Int_t idir = 1;
01860 
01861   ArrivalTime arrtime;
01862 
01863   // determine directionality using timing - this should really
01864   //only be needed for cosmic ray events, as beam events always
01865   //increase in z with increasing time.
01866   
01867   //we could try to determine the timing for each side (E, W) of each
01868   //electronics crate.  that would be the most accurate way, but it 
01869   //would also take a lot of time and iterations.  instead, we can 
01870   //probably get away with just doing a global fit to all datapoints
01871   //since we do have timing constants that take care of the side to side
01872   //and crate to crate offsets.  even if a bit of timing hardware is 
01873   //exchanged, this wont really affect our ability to go back and figure
01874   //out the right calibrations later.
01875   
01876   //make some arrays to do the fit. allow a maximum of 1000 points
01877   //to be included in the fit.  there may be events with more than
01878   //1000 digits, but that should be enough to tell you which way the
01879   //event is going.
01880   Double_t zpos[1000];
01881   Double_t time[1000];
01882   Double_t weight[1000];
01883   
01884   for (int i=0; i<1000; i++) {
01885     zpos[i] = 0.;
01886     time[i] = 0.;
01887     weight[i] = 0.;
01888   }
01889   
01890   //minimum weight for a signal in the timing fit
01891   // maximum weight corresponds to 10 p.e.; reason for this is to minimize
01892   // weight of showers in which transverse spread can have negative effect
01893   // on timing
01894   const Double_t min_wgt = 0.4;
01895   TrackClusterSR *cluster = 0;
01896 
01897   Int_t clusterCtr = 0;
01898 
01899   clusterItr.ResetFirst();
01900   while( (cluster = clusterItr()) && clusterCtr<1000){
01901     
01902     //get the z position of the plane
01903     zpos[clusterCtr] = (Double_t)cluster->GetZPos();
01904         
01905     time[clusterCtr] = cluster->GetBegTime();
01906                 
01907     weight[clusterCtr] = min(arrtime.Weight(cluster->GetCharge()),min_wgt);
01908     ++clusterCtr;
01909 
01910   }//end loop over clusters to fill timing arrays
01911   
01912   //reset the iterator over all strips
01913   clusterItr.ResetFirst();
01914   
01915   //arrays to hold the slope, intercept, and chi^2 of the 
01916   //least squares fit for the timing
01917   Double_t parm[2],eparm[2];
01918   parm[0]=0;
01919   parm[1]=0;
01920   eparm[0]=0;
01921   eparm[1]=0;
01922   //set the maximum residual to the input parameter+1ns to 
01923   //get ready for the while loop below
01924   Double_t maxresid = fParmMaxTimingResid+1.*Munits::ns;
01925   Double_t resid = 0.;
01926   Int_t imaxresid = 0;
01927   Int_t nremove=0;
01928   Int_t fitflag = 0;
01929 
01930   while(clusterCtr-nremove>4 && maxresid>fParmMaxTimingResid && !fitflag){
01931     fitflag = LinearFit::Weighted(clusterCtr,zpos,time,weight,parm,eparm);
01932     maxresid = 0.;
01933     imaxresid = 0;
01934     for (int i=0; i<clusterCtr; ++i) {
01935       resid = parm[0]+parm[1]*zpos[i]-time[i];
01936       if (weight[i]>0. && TMath::Abs(resid)>maxresid) {
01937         maxresid = TMath::Abs(resid);
01938         imaxresid = i;
01939       }
01940     }
01941     if (maxresid>fParmMaxTimingResid) {
01942       weight[imaxresid] = 0.;
01943       ++nremove;
01944     }
01945   }//end loop to fit the timing values
01946 
01947   MSG("AlgTrackSRList", Msg::kDebug) << "cluster timing slope is " << parm[1] << endl;
01948 
01949   if(parm[1]<0.){
01950     idir=-1;
01951   }
01952 
01953   return idir;
01954 }

Int_t AlgTrackSRList::FindTimingDirection ( CandStripHandleItr &  allStripItr,
Float_t  uTrackSlope,
Float_t  uTrackIntercept,
Float_t  vTrackSlope,
Float_t  vTrackIntercept 
) [private]

Definition at line 1691 of file AlgTrackSRList.cxx.

References digit(), fParmMaxTimingResid, fSimFlag, fVldc, CandDigitHandle::GetCharge(), CandHandle::GetDaughterIterator(), PlexSEIdAltL::GetDemuxVetoFlag(), PlexSEIdAltL::GetEnd(), UgliStripHandle::GetHalfLength(), CandStripHandle::GetNDigit(), PlexSEIdAltL::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), UgliGeomHandle::GetScintPlnHandle(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandStripHandle::GetTime(), UgliPlnHandle::GetZ0(), CandStripHandle::GetZPos(), Msg::kDebug, StripEnd::kNegative, CalDigitType::kPE, PlaneView::kU, PlaneView::kV, Msg::kVerbose, min, MSG, Munits::ns, PropagationVelocity::Velocity(), ArrivalTime::Weight(), and LinearFit::Weighted().

Referenced by RunAlg().

01696 {
01697   UgliGeomHandle ugh(fVldc);
01698 
01699   Int_t idir = 1;
01700 
01701   ArrivalTime arrtime;
01702 
01703   // determine directionality using timing - this should really
01704   //only be needed for cosmic ray events, as beam events always
01705   //increase in z with increasing time.
01706   
01707   //we could try to determine the timing for each side (E, W) of each
01708   //electronics crate.  that would be the most accurate way, but it 
01709   //would also take a lot of time and iterations.  instead, we can 
01710   //probably get away with just doing a global fit to all datapoints
01711   //since we do have timing constants that take care of the side to side
01712   //and crate to crate offsets.  even if a bit of timing hardware is 
01713   //exchanged, this wont really affect our ability to go back and figure
01714   //out the right calibrations later.
01715   
01716   //make some arrays to do the fit. allow a maximum of 1000 points
01717   //to be included in the fit.  there may be events with more than
01718   //1000 digits, but that should be enough to tell you which way the
01719   //event is going.
01720   Double_t zpos[1000];
01721   Double_t time[1000];
01722   Double_t weight[1000];
01723   
01724   for (int i=0; i<1000; i++) {
01725     zpos[i] = 0.;
01726     time[i] = 0.;
01727     weight[i] = 0.;
01728   }
01729   
01730   //declare variables used in the while loop below
01731   //distFromCenter is the offset in u or v from the center of the strip
01732   //halflength is the strip halflength
01733   //upos and vpos are obvious
01734   //fiberDist is the total distance traveled in fiber 
01735   Int_t digitCtr = 0;
01736   Double_t distFromCenter = 0.;
01737   Double_t halflength = 0.;
01738   Double_t upos = 0.;
01739   Double_t vpos = 0.;
01740   Double_t fiberDist = 0.;
01741   CandStripHandle *strip = 0;
01742   CandDigitHandle *digit = 0;
01743 
01744   PlaneView::PlaneView_t view0 = PlaneView::kU;
01745   PlaneView::PlaneView_t view1 = PlaneView::kV;
01746 
01747   //minimum weight for a signal in the timing fit
01748   // maximum weight corresponds to 10 p.e.; reason for this is to minimize
01749   // weight of showers in which transverse spread can have negative effect
01750   // on timing
01751   const Double_t min_wgt = 0.4;
01752   
01753   allStripItr.ResetFirst();
01754   MSG("AlgTrackSRList", Msg::kDebug) << "-2 fit over " << digitCtr << " points" << endl;
01755 
01756   while( (strip = allStripItr()) && digitCtr<1000){
01757     
01758     //get the Plex and Ugli information for this strip to be used when 
01759     //figuring out the fiber lengths
01760     PlexStripEndId stripid = strip->GetStripEndId();
01761     UgliScintPlnHandle planehandle = ugh.GetScintPlnHandle(stripid);
01762     TIter digitItr(strip->GetDaughterIterator());
01763     UgliStripHandle striphandle = ugh.GetStripHandle(stripid);
01764     halflength = striphandle.GetHalfLength();
01765     upos = uTrackIntercept+uTrackSlope*strip->GetZPos();
01766     vpos = vTrackIntercept+vTrackSlope*strip->GetZPos();
01767     
01768     //loop over the digits associated with this strip and only use those
01769     //which dont have the demux veto flag set - those that do were not
01770     //used in the demuxing solution
01771     for (int j=0; j<strip->GetNDigit() && digitCtr<1000; ++j) {
01772       digit = dynamic_cast<CandDigitHandle*>(digitItr());
01773       
01774 //------------>>>>coud also require a minimum pulseheight to be included in the fit
01775       if (!digit->GetPlexSEIdAltL().GetDemuxVetoFlag()) {
01776         //get the z position of the plane
01777         zpos[digitCtr] = (Double_t)planehandle.GetZ0();
01778         
01779         PlexSEIdAltL altl = digit->GetPlexSEIdAltL();
01780         StripEnd::StripEnd_t stripEnd = altl.GetEnd();
01781         distFromCenter = 0.;
01782         
01783         // if we're in a U plane we need to account for the V offset
01784         // and vice versa
01785         if (digit->GetPlexSEIdAltL().GetPlaneView() == view1) {
01786           distFromCenter = 
01787             (digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kNegative) ?  upos : -upos;
01788         }
01789         else if (digit->GetPlexSEIdAltL().GetPlaneView() == view0) {
01790           distFromCenter = 
01791             (digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kNegative) ? -vpos :  vpos;
01792         }
01793         
01794         fiberDist = ((halflength + distFromCenter) +
01795                      striphandle.ClearFiber(stripEnd) + 
01796                      striphandle.WlsPigtail(stripEnd));
01797         
01798         time[digitCtr] = (strip->GetTime(digit->GetPlexSEIdAltL().GetEnd()) 
01799                           - (fiberDist/PropagationVelocity::Velocity(fSimFlag)));
01800         
01801         weight[digitCtr] = min(arrtime.Weight(digit->GetCharge(CalDigitType::kPE)),min_wgt);
01802         ++digitCtr;
01803         
01804         MSG("AlgTrackSRList", Msg::kDebug) << "-1 fit over " << digitCtr << " points" << endl;
01805 
01806       }//end if this is a non-vetoed digit
01807     }//end loop over digits
01808   }//end loop over strips to fill timing arrays
01809   
01810   //reset the iterator over all strips
01811   allStripItr.ResetFirst();
01812   
01813   //arrays to hold the slope, intercept, and chi^2 of the 
01814   //least squares fit for the timing
01815   Double_t parm[2],eparm[2];
01816   
01817   //set the maximum residual to the input parameter+1ns to 
01818   //get ready for the while loop below
01819   Double_t maxresid = fParmMaxTimingResid+1.*Munits::ns;
01820   Double_t resid = 0.;
01821   Int_t imaxresid = 0;
01822   Int_t nremove=0;
01823   Int_t fitflag=0;
01824   
01825   MSG("AlgTrackSRList", Msg::kDebug) << "0 fit over " << digitCtr << " points" << endl;
01826   
01827   while(digitCtr-nremove>4 && maxresid>fParmMaxTimingResid && !fitflag){
01828     MSG("AlgTrackSRList", Msg::kDebug) << "1 fit over " << digitCtr << " points" << endl;
01829     fitflag = LinearFit::Weighted(digitCtr,zpos,time,weight,parm,eparm);
01830     maxresid = 0.;
01831     imaxresid = 0;
01832     for (int i=0; i<digitCtr; ++i) {
01833       resid = parm[0]+parm[1]*zpos[i]-time[i];
01834       if (weight[i]>0. && TMath::Abs(resid)>maxresid) {
01835         maxresid = TMath::Abs(resid);
01836         imaxresid = i;
01837       }
01838     }
01839     if (maxresid>fParmMaxTimingResid) {
01840       weight[imaxresid] = 0.;
01841       ++nremove;
01842     }
01843     MSG("AlgTrackSRList", Msg::kDebug) << "2 fit over " << digitCtr << " points" << endl;
01844   }//end loop to fit the timing values
01845 
01846   MSG("AlgTrackSRList", Msg::kVerbose) << "timing slope is " << parm[1] << endl;
01847 
01848   if(parm[1]<0.){
01849     idir=-1;
01850   }
01851 
01852   return idir;
01853 }

void AlgTrackSRList::FormCandTrackSR ( Track2DSRItr &  uTrackItr,
Track2DSRItr &  vTrackItr,
TObjArray *  candTracks,
CandHandle ch,
CandSliceHandle slice,
AlgHandle  ah,
CandContext  cxx,
Int_t  direction 
) [private]

Definition at line 3756 of file AlgTrackSRList.cxx.

References CandHandle::AddDaughterLink(), fDetector, fParmDiffViewBegPlaneMatch, fParmDiffViewEndPlaneMatch, fParmDiffViewPlaneOverlap, fParmDiffViewTimeMatch, fVldc, Track2DSR::GetBegPlane(), Track2DSR::GetCluster(), CandHandle::GetDaughterIterator(), Track2DSR::GetEndPlane(), UgliStripHandle::GetHalfLength(), Track2DSR::GetLast(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), Track2DSR::GetT0(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), UgliStripHandle::GlobalPos(), Msg::kDebug, Detector::kNear, kU, PlaneView::kU, kV, PlaneView::kV, PlaneView::kX, PlaneView::kY, Munits::m, CandTrackSR::MakeCandidate(), MSG, CandRecoHandle::SetCandSlice(), CandContext::SetDataIn(), and tc.

Referenced by RunAlg().

03761 {
03762   Track2DSR *tracku = 0;
03763   Track2DSR *trackv = 0;
03764 
03765   uTrackItr.ResetFirst();
03766   vTrackItr.ResetFirst();
03767   Int_t nUtracks=0;
03768   Int_t nVtracks=0;
03769   while( (tracku = uTrackItr()) )nUtracks++;  
03770   while( (trackv = vTrackItr()) )nVtracks++;
03771   uTrackItr.ResetFirst();
03772   vTrackItr.ResetFirst();
03773 
03774   TrackClusterSR *tc = 0;
03775 
03776   Double_t nhit0;
03777   Bool_t match = false;
03778   Int_t nhitoverlap = 0;
03779   Int_t uTrackCtr = 0, vTrackCtr = 0;
03780 
03781   //loop over the tracks from the 2 views
03782 
03783 
03784   while( (tracku = uTrackItr()) ){
03785 
03786    MSG("TrackSR",Msg::kDebug) << "forming 3D tracks, u track " 
03787                                      << uTrackCtr << " beg plane = " 
03788                                      << tracku->GetBegPlane()
03789                                      << " end plane = " << tracku->GetEndPlane() 
03790                                      << endl;
03791 
03792 
03793    vTrackCtr = 0;
03794 
03795     while( (trackv = vTrackItr()) ){
03796 
03797       match=false;
03798       Int_t begMisses=0;
03799       Int_t endMisses=0;
03800       for(Int_t iu=0; iu<=tracku->GetLast(); ++iu){
03801 
03802         tc = tracku->GetCluster(iu);
03803         if(direction*tc->GetPlane()<direction*trackv->GetBegPlane())begMisses++;
03804         if(direction*tc->GetPlane()>direction*trackv->GetEndPlane())endMisses++; 
03805       }//end loop over u clusters
03806       
03807       for(Int_t iv=0; iv<=trackv->GetLast(); ++iv){
03808         tc = trackv->GetCluster(iv);
03809         if(direction*tc->GetPlane()<direction*tracku->GetBegPlane())begMisses++;
03810         if(direction*tc->GetPlane()>direction*tracku->GetEndPlane())endMisses++; 
03811       }//end loop over v clusters
03812       
03813       
03814       MSG("TrackSR",Msg::kDebug) << " # skip Beg/End " << begMisses << " " << endMisses << endl;
03815 
03816       MSG("TrackSR",Msg::kDebug)  << " u track " << uTrackCtr << "/" << nUtracks  << " v track " << vTrackCtr << "/" << nVtracks  
03817                                         << " beg plane " << trackv->GetBegPlane() 
03818                                        << " end plane " << trackv->GetEndPlane() 
03819                                        << endl
03820                                        << " dtime "
03821                                        << TMath::Abs(tracku->GetT0()
03822                                                      -trackv->GetT0()) 
03823                                        << "/" 
03824                                        << fParmDiffViewTimeMatch << " beg Misses" 
03825                                        << begMisses
03826                                        << "/" << fParmDiffViewBegPlaneMatch 
03827                                        << " end Misses  " 
03828                                        << endMisses 
03829                                        << "/" << fParmDiffViewEndPlaneMatch <<endl;
03830       
03831       // check 2D tracks for position and time agreement. If there
03832       // is only 1 2D track in each view, skip time agreement check
03833       if((TMath::Abs(tracku->GetT0()-trackv->GetT0())
03834           <=fParmDiffViewTimeMatch || 
03835           (nUtracks==1 && nVtracks==1))  
03836          &&(begMisses<=fParmDiffViewBegPlaneMatch 
03837             ||  endMisses <= fParmDiffViewEndPlaneMatch)
03838          && (TMath::Abs(tracku->GetBegPlane()-trackv->GetBegPlane())<= fParmDiffViewBegPlaneMatch*6 || 
03839              TMath::Abs(tracku->GetEndPlane()-trackv->GetEndPlane())<= fParmDiffViewEndPlaneMatch*6)
03840          ){
03841                  
03842         // check hit overlap
03843         nhitoverlap = 0;
03844         //set nhit0 to the minimum of the number of clusters (ie planes) 
03845         //from each view times the amount of overlap required
03846         nhit0 = (1.*TMath::Min(tracku->GetLast()+1,trackv->GetLast()+1)
03847                  *fParmDiffViewPlaneOverlap);
03848         //loop over the clusters in the u track. count up the number of 
03849         //clusters in the u view that are between the first and last planes 
03850         //in the v view track
03851         for(Int_t iu=0; iu<=tracku->GetLast() && 1.*nhitoverlap<nhit0; ++iu){
03852           tc = tracku->GetCluster(iu);
03853           if(direction*tc->GetPlane()>=direction*trackv->GetBegPlane()
03854              && direction*tc->GetPlane()<=direction*trackv->GetEndPlane()) 
03855             ++nhitoverlap;
03856         }//end loop over u clusters
03857         
03858         if(nhitoverlap>=nhit0) match = true;
03859 
03860         if(!match){
03861           nhitoverlap = 0;
03862           //no match found, loop over v clusters
03863           for(Int_t iv=0; iv<=trackv->GetLast() && 1.*nhitoverlap<nhit0; ++iv){
03864             tc = trackv->GetCluster(iv);
03865             if(direction*tc->GetPlane()>=direction*tracku->GetBegPlane()
03866                && direction*tc->GetPlane()<=direction*tracku->GetEndPlane()) 
03867               ++nhitoverlap;
03868           }//end loop over v clusters
03869 
03870           if (nhitoverlap>=nhit0) match = true;
03871 
03872         }//end if !match
03873        
03874 
03875         //does it match now?, if so, make a CandTrackSRHandle
03876         if(match){
03877           TObjArray trackin;
03878           trackin.Add(tracku);
03879           trackin.Add(trackv);
03880           trackin.Add(slice);
03881           cxx.SetDataIn(&trackin);
03882           CandTrackSRHandle trackhandle = CandTrackSR::MakeCandidate(ah,cxx);
03883           trackhandle.SetCandSlice(slice);
03884           // count number of track hits with 3D location outside detector
03885           UgliGeomHandle ugh(fVldc);
03886           CandStripHandleItr stripItr(trackhandle.GetDaughterIterator());
03887           Float_t nstrip=0;
03888           Float_t noutside=0;
03889           CandStripHandle *strip=0;
03890           while( (strip = stripItr()) ){
03891             PlexStripEndId stripid = strip->GetStripEndId();
03892             UgliStripHandle striphandle = ugh.GetStripHandle(stripid);
03893             Double_t halflength = striphandle.GetHalfLength();
03894             Int_t plane=strip->GetPlane();
03895             Double_t stripX=0;
03896             switch (strip->GetPlaneView()) {
03897             case PlaneView::kX: case PlaneView::kU:
03898               stripX = - trackhandle.GetV(plane);
03899               break;
03900             case PlaneView::kY: case PlaneView::kV:
03901               stripX = trackhandle.GetU(plane);
03902               break;
03903             default:
03904               continue;
03905             }
03906             if(fDetector==Detector::kNear) {
03907               const double rsqrt2 = 1.0/sqrt(2.0);
03908               const TVector3 kU(rsqrt2,rsqrt2,0.);
03909               const TVector3 kV(-rsqrt2,rsqrt2,0.);  
03910               if (strip->GetPlaneView()==PlaneView::kU) {
03911                 double lambda = kV.Dot(striphandle.GlobalPos(0.));
03912                 stripX += lambda;
03913               } else {
03914                 double lambda = kU.Dot(striphandle.GlobalPos(0.));
03915                 stripX -= lambda;
03916               }
03917             }
03918             if(TMath::Abs(stripX)>halflength+1.0*Munits::m){
03919               noutside++;
03920             }
03921             nstrip++;
03922           }
03923           
03924           // if less than 3/4 of hit planes are within active detector volume, don't keep this track
03925            if(noutside/nstrip<0.25){
03926             // the following should be cast to CandTrackSRHandle, but doing 
03927             // so gives a compilation error
03928             const CandTrackSRHandle *track = dynamic_cast<const CandTrackSRHandle*>
03929               (ch.AddDaughterLink(trackhandle));
03930             candTracks->AddLast(const_cast<CandTrackSRHandle*>(track));
03931               }
03932            else{
03933              //      delete trackhandle;
03934              MSG("TrackSR", Msg::kDebug) << " **** removing track with " << noutside/nstrip << " outside detector " << endl;
03935            }
03936          
03937         }
03938 //      else {
03939 //        MSG("TrackSR",Msg::kDebug) << "no match, plane overlap = " 
03940 //                                   << nhitoverlap 
03941 //                                   << "/" << nhit0 << "\n";
03942 //      }
03943       }//end if tracks could match based on extent
03944       ++vTrackCtr;
03945     }//end loop over v tracks
03946    
03947     //reset the v track iterator
03948     vTrackItr.ResetFirst();
03949 
03950     ++uTrackCtr;
03951   }//end loop over u tracks
03952 
03953   MSG("TrackSR", Msg::kDebug) << candTracks->GetLast()+1 
03954                               << " tracks formed" << endl;
03955 
03956   return;
03957 }

Int_t AlgTrackSRList::IdentifyBadTracks ( Track2DSRItr &  trackItr,
Int_t  iterate,
Int_t  trk2dmin 
) [private]

Definition at line 2921 of file AlgTrackSRList.cxx.

References Track2DSR::fIterate, fParmTrk2DHitFraction, fParmTrk2DMaxResid, fSliceDzDs, Track2DSR::GetBegPlane(), Track2DSR::GetCluster(), Track2DSR::GetEndPlane(), Track2DSR::GetLast(), TrackClusterSR::GetPlane(), TrackClusterSR::GetTPos(), TrackClusterSR::GetZPos(), Track2DSR::IsBad(), Msg::kDebug, max, MSG, Track2DSR::RemoveAt(), Track2DSR::SetIsBad(), and LinearFit::Unweighted().

Referenced by RunAlg().

02923 {
02924   //get to the first track in the list
02925   trackItr.ResetFirst();
02926   Track2DSR *track = 0;
02927   TrackClusterSR *tctmp = 0;
02928   Double_t hitfrac = 0.;
02929   Int_t ilast = 0;
02930   Int_t nhit = 0;
02931   Int_t removedTracks = 0; // Double_t dplane = 0.;
02932   while( (track = trackItr()) ){ //loop over the tracks
02933 
02934     //if the track was created in this iteration and isnt already bad
02935     if(track->fIterate==iterate && !track->IsBad() ){
02936       hitfrac = 0.;
02937       Bool_t hitRemoved=true;
02938       while(hitRemoved && track->GetLast()>0 ){
02939         ilast = track->GetLast();
02940         nhit = ilast+1;
02941         hitRemoved=false;
02942         Float_t aveSpacing= 0.5*TMath::Abs(track->GetEndPlane()-
02943                                    track->GetBegPlane())/ilast;
02944         Float_t lastSpacing=0.5*TMath::Abs(track->GetCluster(ilast)->GetPlane()-
02945                                              track->GetCluster(ilast-1)->GetPlane());
02946         hitfrac=aveSpacing/lastSpacing;
02947 
02948         
02949         //remove the last clusters from tracks with too few hits for the number 
02950         //of planes crossed until the hit occupancy is acceptable
02951         if(hitfrac<fParmTrk2DHitFraction && lastSpacing>5){
02952           track->RemoveAt(ilast);
02953           hitRemoved=true;
02954           MSG("TrackSR",Msg::kDebug) << "removing end hit: hit fraction " 
02955                                      << hitfrac << "/" << fParmTrk2DHitFraction 
02956                                      <<endl;
02957         }
02958 
02959         
02960       }//end loop while not enough hit occupancy
02961         
02962       //remove short tracks
02963       if(track->GetLast()+1<trk2dmin){
02964         ++removedTracks;
02965         track->SetIsBad(true);
02966         MSG("TrackSR",Msg::kDebug) << "adding badtrack: short track " 
02967                                    << track->GetLast()+1 << "/" << trk2dmin 
02968                                    << " beg plane = " << track->GetBegPlane()
02969                                    << endl;
02970       }//end if short track
02971 
02972 
02973       //check linearity of short tracks
02974       nhit = track->GetLast()+1;
02975       if(nhit<=fParmTrk2DNHough && nhit>2){
02976         Double_t *xfit = new Double_t[nhit];
02977         Double_t *yfit = new Double_t[nhit];
02978         Double_t parm[2];
02979         for (Int_t ifit=0; ifit<nhit; ++ifit) {
02980           tctmp = track->GetCluster(ifit);
02981           xfit[ifit] = tctmp->GetZPos();
02982           yfit[ifit] = tctmp->GetTPos();
02983         }
02984         LinearFit::Unweighted(nhit,xfit,yfit,parm);
02985         Double_t maxresid = 0.;
02986         for (Int_t ifit=0; ifit<nhit; ++ifit) {
02987           maxresid = max(maxresid, TMath::Abs(parm[0]+parm[1]*xfit[ifit]-yfit[ifit]));
02988         }
02989         delete [] xfit;
02990         delete [] yfit;
02991         if(maxresid>fParmTrk2DMaxResid*max(fSliceDzDs,
02992                                            TMath::Sqrt(1.+parm[1]*parm[1]))){
02993           ++removedTracks;
02994           track->SetIsBad(true);
02995           MSG("TrackSR",Msg::kDebug) << "adding badtrack: linearity " 
02996                                      << maxresid << "/" 
02997                                      << fParmTrk2DMaxResid <<endl;
02998         }
02999 
03000       }//end if short track
03001 
03002     }//end if track created in current iteration
03003 
03004     MSG("TrackSR", Msg::kDebug) << "track with beg plane " << track->GetBegPlane()
03005                                << " is bad " << (Int_t)track->IsBad() << endl;
03006   }//end loop over tracks
03007  
03008   return removedTracks;
03009 }

Bool_t AlgTrackSRList::IsBestClusterInPlane ( TrackClusterSR planeCluster,
TrackClusterSR prevCluster,
Int_t  numSkippedHit,
Int_t  nHit,
Int_t  planesBefore,
Double_t  trackSlope,
Track2DSR track,
Double_t *  residParams,
Int_t  dplane,
Bool_t  isSpectrometerPlane 
) [private]

Definition at line 2616 of file AlgTrackSRList.cxx.

References Mphysical::c_light, fDetector, fParmHitNTime, fParmHitTime, fParmIsCosmic, fParmMinSingleHit, fParmSingleHitDef, fParmTrk2DDirCosNSigma, fParmTrk2DLinA0, fParmTrk2DLinB0, fParmTrk2DNHough, fParmTrk2DNHough0, fParmTrk2DNSigmaA, fParmTrk2DNSkipHit, fParmTrkNPlaneGoodDir, fSliceEndPlane, fSliceNUPlanes, fSliceNVPlanes, fSliceVtxPlane, TrackClusterSR::GetBegTime(), TrackClusterSR::GetCharge(), Track2DSR::GetCluster(), Track2DSR::GetHoughExist(), Track2DSR::GetHoughSlope(), TrackClusterSR::GetMaxTPos(), TrackClusterSR::GetMinTPos(), TrackClusterSR::GetNStrip(), TrackClusterSR::GetPlane(), Track2DSR::GetT0(), TrackClusterSR::GetTPos(), TrackClusterSR::GetTPosError(), Track2DSR::GetZ0(), TrackClusterSR::GetZPos(), Msg::kDebug, Detector::kFar, Detector::kNear, Munits::MeV, min, MSG, and LinearFit::Weighted().

Referenced by AddClustersToTracks(), and SpectrometerTracking().

02624 {
02625   
02626   bool bestCluster = false;
02627   bool useph = (fDetector == Detector::kFar && fParmIsCosmic==1);
02628 
02629  // various constants used in this method
02630   Double_t sigmams = 0.014;
02631   Double_t radLenPerPlane = 1.4;
02632   Double_t dPPerPlane=40*Munits::MeV;
02633   Double_t nomBField=0.8; // Telsa
02634 
02635   //get a pointer to a cluster
02636   TrackClusterSR *tctmp = 0;
02637 
02638   //get variables to describe the differences from the prevCluster to the 
02639   //one in the current plane in z, time, transverse position,etc
02640   Double_t dzpos = TMath::Abs(planeCluster->GetZPos()-prevCluster->GetZPos());
02641   Double_t dtime = planeCluster->GetBegTime()-track->GetT0()-
02642     (planeCluster->GetZPos()-track->GetZ0())/Mphysical::c_light;
02643    
02644   Double_t resid=0.;
02645   Double_t maxresid=99999.;
02646   Double_t dangle = 0.;
02647   Double_t maxdangle = 999999.;
02648 
02649   //declare the size of the window in # of planes over which you want to 
02650   //find the slope of the track
02651   // init maximum allowable residual fixed sized (related to strip width) plus
02652   // contribution scaling with distance extrapolated
02653   Int_t mhit;
02654   if(nHit<fParmTrk2DNHough0){
02655     mhit = nHit;    
02656     maxresid = fParmTrk2DLinA0+fParmTrk2DLinB0*TMath::Abs(dzpos);
02657   }
02658   else{
02659     mhit = min(nHit,fParmTrk2DNHough);
02660     if( fDetector == Detector::kNear && prevCluster->GetPlane()>121 && mhit>2)mhit--;
02661     maxresid = TMath::Abs(planeCluster->GetMaxTPos()-planeCluster->GetMinTPos())*0.288;
02662   }
02663   maxresid *=maxresid;
02664   if(nHit>1){     
02665 
02666     //nexitplane tells you how many planes you have left between the 
02667     //plane of the last cluster added to the track and the most downstream hit in the
02668     //slice. This is a (crude) estimate of the muon energy at that point,  and
02669     // allows an estimate of bending and MCS 
02670     Int_t nexitplane = TMath::Abs(prevCluster->GetPlane()-fSliceEndPlane)+1;
02671     if(fSliceNUPlanes+fSliceNVPlanes<fParmTrkNPlaneGoodDir && fParmIsCosmic==1){
02672       // this is a short track, so the timing is not so good
02673       // take smaller momentum from either end
02674       nexitplane = min(nexitplane,TMath::Abs(prevCluster->GetPlane()-fSliceVtxPlane)+1);
02675     }
02676        
02677     Double_t *xfit = new Double_t[mhit];
02678     Double_t *yfit = new Double_t[mhit];
02679     Double_t *wfit = new Double_t[mhit];
02680     Double_t parm[2];
02681     Double_t eparm[2] = {0.,0.};
02682     Double_t wtmp = 0.;
02683  
02684     //fit a straight line over the clusters in a window of size
02685     //mhit to see how well the current cluster matches up to the
02686     //fit of previously added clusters
02687     for (Int_t ifit=0; ifit<mhit; ++ifit) {
02688       if(isSpectrometerPlane){
02689         tctmp = track->GetCluster(ifit+nBef);
02690       }
02691       else{
02692         tctmp = track->GetCluster(nHit-mhit+ifit);
02693       }
02694       xfit[ifit] = tctmp->GetZPos()-planeCluster->GetZPos();
02695       yfit[ifit] = tctmp->GetTPos();
02696       wtmp = 0.288*(tctmp->GetMaxTPos()-tctmp->GetMinTPos());
02697       wfit[ifit]=1.0/(wtmp*wtmp);
02698     }
02699 
02700     for(Int_t ifit=0;ifit<mhit;ifit++){
02701       MSG("TrackSR", Msg::kDebug)  << " fit input " << xfit[ifit]  << " " << yfit[ifit] << " " << wfit[ifit] << endl; 
02702     }
02703     
02704     LinearFit::Weighted(mhit,xfit,yfit,wfit,parm,eparm);
02705     
02706     delete [] xfit;
02707     delete [] yfit;
02708     delete [] wfit;
02709 
02710     //get the residual for this plane to the fit over the previous planes
02711     MSG("TrackSR", Msg::kDebug) << "fit slope = " << parm[1] <<  "/ " << eparm[1] 
02712                                 << " fit intercept = " 
02713                                 << parm[0] << " / " << eparm[0] << endl;
02714 
02715     //find the change in the angle from the fit slope 
02716     Double_t planeClusterSlope = (planeCluster->GetTPos()-prevCluster->GetTPos())
02717                                   /(planeCluster->GetZPos()-prevCluster->GetZPos());
02718 
02719     dangle = TMath::Abs(parm[1]-planeClusterSlope);
02720 
02721     resid = TMath::Abs(parm[0]-planeCluster->GetTPos());
02722 
02723     // determine max allowable residuals
02724     //get the range-based  momentum assuming a loss of 40 MeV per plane
02725     Double_t dzds = 1.0/TMath::Sqrt(1.+trackSlope*trackSlope);
02726     Double_t pathlength = (Double_t)(nexitplane)/dzds;
02727     Double_t tctp = pathlength*dPPerPlane;
02728     Double_t dbfangle =  0.3*dzpos/dzds*nomBField/tctp;
02729     Double_t mcsangle =  sigmams*TMath::Sqrt(radLenPerPlane*TMath::Abs(dplane)/dzds)/tctp;     
02730 
02731     maxdangle = planeCluster->GetTPosError()*planeCluster->GetTPosError()/(dzpos*dzpos) +
02732       dbfangle*dbfangle + mcsangle*mcsangle + eparm[1]*eparm[1] ; 
02733     maxresid += dzpos*dzpos*(dbfangle*dbfangle + mcsangle*mcsangle) +  eparm[0]*eparm[0]; 
02734  
02735     maxresid  =  fParmTrk2DNSigmaA*TMath::Sqrt(maxresid);
02736     maxdangle =  fParmTrk2DDirCosNSigma*TMath::Sqrt(maxdangle);
02737 
02738     // make requirements more stringent if skipping already found clusters
02739     if(nBef>0){
02740       maxresid /=2.;
02741       maxdangle /=2.;
02742     }
02743  
02744     MSG("TrackSR",Msg::kDebug) << " max resids:  dzpos " << dzpos 
02745                                << " bfield " << dbfangle 
02746                                << " mcs " << mcsangle 
02747                                << " maxdangle " << maxdangle 
02748                                << " bfield resid " << dbfangle*dzpos 
02749                                << " mcs resid " << mcsangle*dzpos << endl;    
02750   }//end if at least 2 hits
02751   
02752   else if(track->GetHoughExist()){
02753     //end not enough hits, try getting residuals using the hough transform values 
02754     tctmp = track->GetCluster(0);
02755     resid = min(TMath::Abs((tctmp->GetTPos()+(planeCluster->GetZPos()-tctmp->GetZPos())*track->GetHoughSlope())-planeCluster->GetTPos()),
02756                 TMath::Abs(tctmp->GetTPos()-planeCluster->GetTPos()));
02757    
02758     maxresid  =  fParmTrk2DNSigmaA*TMath::Sqrt(maxresid +
02759                                                (tctmp->GetMaxTPos()-tctmp->GetMinTPos())*
02760                                                (tctmp->GetMaxTPos()-tctmp->GetMinTPos())*0.083);  
02761   }
02762 
02763   //get the weight residual based on an empirical parameterization done by
02764   //R. Lee
02765   Double_t weightedresid = (0.3+exp(-0.2*planeCluster->GetCharge()))*resid;
02766   maxresid = min(0.4,maxresid);  // don't allow maxresid to exceed 0.4 m
02767   
02768   MSG("TrackSR",Msg::kDebug) << "plane " << planeCluster->GetPlane() << "numskiphit " <<
02769     numSkippedHit << "/" << fParmTrk2DNSkipHit << endl <<
02770     " nhit " << nHit << "/" << fParmMinSingleHit << " nstrip " << planeCluster->GetNStrip() << "/" << fParmSingleHitDef <<  endl <<
02771     " dtime " << dtime*1.e9 << "/" << fParmHitNTime*1.e9 << "/" << fParmHitTime*1.e9 <<  endl << 
02772     " resid " << resid << "/" << maxresid  << " dangle " << dangle << "/" << maxdangle << endl <<
02773     " useph " << useph << "/" << residParams[0]  << "/" << TMath::Abs(weightedresid) << "/" << residParams[1] << endl; 
02774 
02775 
02776   // criteria for this cluster being the best cluster on this plane to add to track
02777   if(numSkippedHit<=fParmTrk2DNSkipHit 
02778      && (nHit>fParmMinSingleHit || planeCluster->GetNStrip()<=fParmSingleHitDef) 
02779      && TMath::Abs(resid)<maxresid && dangle<maxdangle
02780      && dtime>fParmHitNTime && dtime<fParmHitTime 
02781      && ((!useph && TMath::Abs(resid)<residParams[0]) 
02782          || (useph && TMath::Abs(weightedresid)<residParams[1]))){
02783     MSG("TrackSR",Msg::kDebug) << "  accepted with resid = " <<  resid
02784                                << " weighted resid = "<< weightedresid << endl;
02785     // add hit to this track
02786     bestCluster = true;
02787     residParams[0] = TMath::Abs(resid);
02788     residParams[1] = TMath::Abs(weightedresid);
02789   }//end if this cluster works
02790   
02791   return bestCluster;
02792 }

Int_t AlgTrackSRList::LookForBadClusters ( TrackClusterSRItr &  clusterItr,
TObjArray *  trackList,
Int_t  direction 
) [private]
void AlgTrackSRList::MakeTrackClusters ( TObjArray *  trackClusterList,
CandStripHandleItr &  stripItr,
CandTrackSRListHandle cth,
Int_t  expectedStrips,
Int_t  direction 
) [private]

Definition at line 4274 of file AlgTrackSRList.cxx.

References TrackClusterSR::AddStrip(), CandTrackSRListHandle::AddTrackCluster(), fDetector, fParmMaxNExtraStrip, fParmMinClusterPulseHeight, fParmMisalignmentError, fParmTrkClsNSkip, CandStripHandle::GetBegTime(), TrackClusterSR::GetCharge(), CandStripHandle::GetCharge(), TrackClusterSR::GetNStrip(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), CandStripHandle::GetStrip(), TrackClusterSR::GetTPos(), TrackClusterSR::GetZPos(), TrackClusterSR::InShowerLikePlane(), TrackClusterSR::IsWide(), Msg::kDebug, Detector::kNear, MSG, TrackClusterSR::SetInShowerLikePlane(), TrackClusterSR::SetIsWide(), and tc.

Referenced by RunAlg().

04279 {
04280   //array of flags to mark showerlike planes - those with at least one
04281   //wide cluster
04282   Int_t showerLike[500];
04283   for(Int_t i=0; i<500; ++i){
04284     showerLike[i] = 0;
04285   }
04286 
04287   MSG("AlgTrackSRList", Msg::kDebug) << "expectedStrips " << expectedStrips
04288                                      << " direction " << direction << endl;
04289 
04290   Int_t oldplane=0;
04291   Int_t oldstrip=0;
04292   Double_t oldtime=0.;
04293   TrackClusterSR *tc=0;
04294   stripItr.SetDirection(direction);
04295   stripItr.Reset();
04296 
04297   //loop over the strips in the Iterator
04298   //declare a pointer and zero it so you dont have to keep declaring 
04299   //the strip handle
04300   CandStripHandle *strip = 0;
04301 
04302   while((strip = stripItr())){
04303 
04304     // don't include spectrometer hits
04305     if(fDetector!=Detector::kNear || strip->GetPlane()<121){
04306       MSG("TrackSR", Msg::kDebug) << "looking at plane " << strip->GetPlane()
04307                                   << " strip " << strip->GetStrip() << endl; 
04308       
04309       //check to see if you have a current cluster 
04310       if(tc){
04311         //now see if this strip comes from the same plane and 
04312         //is part of the current cluster (ie not separated by too many strips
04313         if(strip->GetPlane()==oldplane &&
04314            direction*strip->GetStrip()<=direction*oldstrip+1+fParmTrkClsNSkip){
04315           tc->AddStrip(strip);
04316           oldstrip = strip->GetStrip();
04317           oldtime = strip->GetBegTime();
04318           
04319           MSG("TrackSR", Msg::kDebug) << "add strip " << oldstrip 
04320                                       << " with time "  << oldtime*1.e9 
04321                                       << " and charge " << strip->GetCharge() 
04322                                       << " to cluster on plane " 
04323                                       << oldplane << " " << fParmTrkClsNSkip 
04324                                       << endl;
04325           
04326         }
04327         else{
04328           //you have a new plane or cluster - see if the current
04329           //cluster passes the pusle height cut before adding it to the array
04330           MSG("TrackSR", Msg::kDebug) << "checking tc charge " << tc->GetCharge() << "/" << fParmMinClusterPulseHeight << endl;
04331         
04332 
04333           if(tc->GetCharge()>=fParmMinClusterPulseHeight){
04334             
04335             //check to see if this is a wide cluster - ie it has more strips
04336             //in it than expected based on the hough transform - the number of 
04337             //expected strips is passed into the method.  if it is wide, set the 
04338             //cluster's wide flag. also set the flag in the showerLike array
04339             //for this plane
04340             if(tc->GetNStrip()>fParmMaxNExtraStrip+expectedStrips){
04341               tc->SetIsWide(true);
04342               showerLike[tc->GetPlane()] = 1;
04343               
04344               MSG("TrackSR", Msg::kDebug) << "setting cluster on plane " 
04345                                           << tc->GetPlane() << " wide with "
04346                                           << tc->GetNStrip() << " strips " 
04347                                           << "expected " << expectedStrips << endl;
04348               
04349             }
04350             
04351             MSG("TrackSR", Msg::kDebug) << "adding new cluster on plane " 
04352                                         << tc->GetPlane() << " to list " 
04353                                         << tc->GetTPos() << " " 
04354                                         << tc->GetZPos() << " with charge " << tc->GetCharge()
04355                                         << endl;
04356             
04357             trackClusterList->AddLast(tc);
04358             cth.AddTrackCluster(tc);
04359           }//end if enough charge in this cluster to add it to the list
04360           else{
04361             delete tc;
04362           }
04363           //now make a new cluster
04364           TrackClusterSR *trackcluster = new TrackClusterSR(strip, 
04365                                                             fParmMisalignmentError);
04366           oldplane = strip->GetPlane();
04367           oldstrip = strip->GetStrip();
04368           oldtime = strip->GetBegTime();
04369           tc = trackcluster;
04370           
04371           MSG("TrackSR", Msg::kDebug) << "make cluster with strip " << oldstrip 
04372                                       << " with time " << oldtime*1e9 
04373                                       << " with charge " << strip->GetCharge() 
04374                                       << " on plane " << oldplane << endl;
04375           
04376         }
04377       }//end if you have a current cluster
04378       else {//make a new cluster
04379         TrackClusterSR *trackcluster = new TrackClusterSR(strip, 
04380                                                           fParmMisalignmentError);
04381         oldplane = strip->GetPlane();
04382         oldstrip = strip->GetStrip();
04383         oldtime = strip->GetBegTime();
04384         tc = trackcluster;
04385         
04386         MSG("TrackSR", Msg::kDebug) << "make cluster with strip " << oldstrip 
04387                                     << " with time " << oldtime*1e9  
04388                                     << " with charge " << strip->GetCharge() 
04389                                     << " on plane " << oldplane << endl;
04390         
04391       }
04392     }
04393   }//end loop over strips
04394 
04395 
04396   //done with loop, but last cluster needs to be added
04397   if(tc){
04398     if(tc->GetCharge()>=fParmMinClusterPulseHeight){
04399       
04400       //check to see if this is a wide cluster - ie it has more strips
04401       //in it than expected based on the hough transform - the number of 
04402       //expected strips is passed into the method.  if it is wide, set the 
04403       //cluster's wide flag. also set the flag in the showerLike array
04404       //for this plane
04405       if(tc->GetNStrip()>fParmMaxNExtraStrip+expectedStrips){
04406         tc->SetIsWide(true);
04407         showerLike[tc->GetPlane()] = 1;
04408       }
04409       
04410       MSG("TrackSR", Msg::kDebug) << "adding new cluster on plane " 
04411                                   << tc->GetPlane() << " to list " 
04412                                   << tc->GetTPos() << " " << tc->GetZPos() << endl;
04413       
04414       trackClusterList->AddLast(tc);
04415       cth.AddTrackCluster(tc);
04416     }
04417     else{
04418       delete tc;
04419     }
04420   }//end if still a cluster to add to the list
04421 
04422   //loop over the clusters in the array to mark ones in showerlike planes
04423   for(Int_t j=0; j<trackClusterList->GetLast(); ++j){
04424     tc = dynamic_cast<TrackClusterSR *>(trackClusterList->At(j)); 
04425     if(showerLike[tc->GetPlane()] == 1) 
04426       tc->SetInShowerLikePlane(true);
04427   }//end loop over clusters 
04428   
04429   //debugging output for the cluster list made
04430    for (Int_t i=0; i<=trackClusterList->GetLast(); ++i) {
04431      tc = dynamic_cast<TrackClusterSR*>(trackClusterList->At(i));
04432      MSG("TrackSR",Msg::kDebug) << "trackcluster " << i << " plane " 
04433                                 << tc->GetPlane() << " tpos " << tc->GetTPos() 
04434                                 << " zpos " << tc->GetZPos() << endl
04435                                 << " nstrip " << tc->GetNStrip() << " charge " 
04436                                 << tc->GetCharge() << " wide " 
04437                                 << (Int_t)tc->IsWide()
04438                                 << " showerlike " 
04439                                 << (Int_t)tc->InShowerLikePlane() 
04440                                 << endl;
04441    }
04442   
04443   
04444   stripItr.ResetFirst();
04445   return;
04446 }

void AlgTrackSRList::MarkUsedClusters ( Track2DSRItr &  trackItr,
TrackClusterSRItr &  clusterItr,
Int_t  iterate,
Int_t  nmaxtrack,
Double_t  houghSlope,
Double_t  houghIntercept 
) [private]

Definition at line 3014 of file AlgTrackSRList.cxx.

References MuELoss::e, Track2DSR::fIterate, TrackClusterSR::GetBegTime(), Track2DSR::GetCluster(), Track2DSR::GetLast(), TrackClusterSR::GetPlane(), TrackClusterSR::GetTPos(), TrackClusterSR::GetZPos(), Track2DSR::IsBad(), TrackClusterSR::IsValid(), Msg::kDebug, MSG, TrackClusterSR::SetIsValid(), and tc.

Referenced by RunAlg().

03018 {
03019   Track2DSR *track = 0;
03020   TrackClusterSR *tc = 0;
03021   TrackClusterSR *tc1 = 0;
03022   TrackClusterSR *tcbeg = 0;
03023   TrackClusterSR *tcend = 0;
03024 
03025   Double_t residbeg =0.;
03026   Double_t residend =0.;
03027 
03028   trackItr.ResetFirst();
03029   clusterItr.ResetFirst();
03030 
03031   //loop over the tracks
03032   while( (track = trackItr()) ){
03033     
03034     //if this track was created in this iteration
03035     if(track->fIterate==iterate){
03036 
03037       //if it was a good track
03038       if( !track->IsBad() ){
03039         MSG("TrackSR",Msg::kDebug) << "removing hits from track\n";
03040 
03041         //loop over the clusters in this track
03042         for(Int_t itc=0; itc<=track->GetLast(); ++itc){
03043           tc = track->GetCluster(itc);
03044           MSG("TrackSR",Msg::kDebug) << "  TC " << tc->GetPlane() << " " 
03045                                        << tc->GetTPos() <<endl;
03046 
03047           //loop over all clusters in the clusterItr and mark all the ones
03048           //used in this track as non-valid
03049           while( (tc1 = clusterItr()) ){
03050 
03051             //if tc1 is valid and the same cluster as tc, mark it as non-valid
03052             //because it has been used
03053             if(tc1->IsValid() && tc1->GetPlane()==tc->GetPlane() 
03054                && TMath::Abs(tc1->GetTPos()-tc->GetTPos())<1.e-3
03055                && TMath::Abs(1.e9*(tc1->GetBegTime()-tc->GetBegTime()))<1.e-3){
03056               
03057               tc1->SetIsValid(false);
03058 
03059               MSG("TrackSR",Msg::kDebug) << "    REMOVING\n";
03060             }//end if the same cluster and valid
03061 
03062           }//end loop over clusters
03063 
03064           clusterItr.ResetFirst();
03065         }//end loop over clusters in this track
03066       }//end if good track
03067       else{ //now do the bad tracks
03068 
03069         if(nmaxtrack>0){
03070           tcbeg = track->GetCluster(0);
03071           tcend = track->GetCluster(track->GetLast());
03072      
03073           residbeg = TMath::Abs(houghIntercept
03074                                 +houghSlope*tcbeg->GetZPos()-tcbeg->GetTPos());
03075           residend = TMath::Abs(houghIntercept
03076                                 +houghSlope*tcend->GetZPos()-tcend->GetTPos());
03077           MSG("TrackSR",Msg::kDebug) << "bad tracks, hough track exists: tcbeg " 
03078                                      << tcbeg->GetPlane() << " " 
03079                                      << tcbeg->GetTPos() 
03080                                      << " " << residbeg << "   tcend " 
03081                                      << tcend->GetPlane() << " " 
03082                                      << tcend->GetTPos() 
03083                                      << " " << residend << endl;
03084 
03085           //remove the cluster with the larger residual
03086           if(residbeg<residend){
03087             MSG("TrackSR",Msg::kDebug) << "  removing end cluster\n";
03088             tc = tcend;
03089 
03090             //loop over all clusters in the clusterItr and mark all the ones
03091             //used in this track as non-valid
03092             while( (tc1 = clusterItr()) ){
03093             
03094               if (tc1->IsValid() && tc1->GetPlane()==tc->GetPlane() &&
03095                   TMath::Abs(tc1->GetTPos()-tc->GetTPos())<1.e-3 &&
03096                   TMath::Abs(1.e9*(tc1->GetBegTime()-tc->GetBegTime()))<1.e-3) {
03097                 tc1->SetIsValid(false);
03098                 MSG("TrackSR",Msg::kDebug) << "  REMOVING HIT\n";
03099               }//end if the same cluster
03100             }//end loop over clusters
03101             clusterItr.ResetFirst();
03102             
03103           } //end if removing end cluster
03104           else{ //removing begining cluster
03105             MSG("TrackSR",Msg::kDebug) << "  removing beg cluster\n";
03106             tc = tcbeg;
03107           }
03108         }//end if at least one track in this view
03109         else{
03110           tc = track->GetCluster(0);
03111           MSG("TrackSR",Msg::kDebug) << "bad tracks, no hough track: tcbeg " 
03112                                        << tc->GetPlane() << " " << tc->GetTPos() 
03113                                        << " " 
03114                                        << tc->GetBegTime()*1.e9 << "\n";
03115         }
03116         
03117         //loop over all clusters in the clusterItr and mark all the ones
03118         //used in this track as non-valid
03119         while( (tc1 = clusterItr()) ){
03120           if (tc1->IsValid() && tc1->GetPlane()==tc->GetPlane() &&
03121               TMath::Abs(tc1->GetTPos()-tc->GetTPos())<1.e-3 &&
03122               TMath::Abs(1.e9*(tc1->GetBegTime()-tc->GetBegTime()))<1.e-3) {
03123             tc1->SetIsValid(false);
03124             MSG("TrackSR",Msg::kDebug) << "  REMOVING HIT\n";
03125           }//end if the same cluster
03126         }//end loop over clusters
03127         clusterItr.ResetFirst();
03128 
03129       }//end else bad track
03130     }//end if track created in this iteration
03131   }//end loop over tracks
03132   
03133   return;
03134 }

void AlgTrackSRList::MergeTracks ( TObjArray *  tracklist,
Int_t  direction 
) [private]

Definition at line 3303 of file AlgTrackSRList.cxx.

References Track2DSR::Add(), Track2DSR::Compress(), fHoughUIntercept, fHoughUSlope, fHoughVIntercept, fHoughVSlope, FindNumSkippedPlanes(), fParmTrk2DDirCosError2, fParmTrk2DLinA, fParmTrk2DLinB, fParmTrk2DNSeed, fParmTrk2DNSkipHit, Track2DSR::GetBegPlane(), Track2DSR::GetCluster(), Track2DSR::GetForwardSlope(), Track2DSR::GetLast(), TrackClusterSR::GetPlane(), VHS::GetPlane(), Track2DSR::GetPlaneView(), TrackClusterSR::GetTPos(), TrackClusterSR::GetTPosError(), TrackClusterSR::GetZPos(), Track2DSR::IsBad(), Msg::kDebug, PlaneView::kU, max, min, MSG, Track2DSR::RemoveAt(), and Track2DSR::SetIsBad().

Referenced by RunAlg().

03305 {
03306   MSG("TrackSR",Msg::kDebug) << " in MergeTracks" << endl;
03307   Track2DSRItr trackItr1(tracklist);  
03308   Track2DSRItr trackItr2(tracklist);
03309   Track2DSR *track1 = 0;
03310   Track2DSR *track2 = 0;
03311   for (Int_t nremove1=0;nremove1<5;nremove1++){ 
03312     for (Int_t nremove2=0;nremove2<5;nremove2++){
03313       Int_t fwd1,fwd2,bck1,bck2;
03314       Bool_t merge=false;
03315       Int_t itrack1=0;
03316       trackItr1.Reset();
03317       while( (track1 = trackItr1()) ){
03318           MSG("TrackSR",Msg::kDebug) <<  " track 1 " << track1->GetCluster(0)->GetPlane() << "/" << 
03319                     track1->GetCluster(track1->GetLast())->GetPlane() << endl;
03320         if(track1->GetLast()>fParmTrk2DNSeed && !track1->IsBad()){
03321           Int_t itrack2=0;
03322           trackItr2.Reset();
03323           while((track2=trackItr2())){ 
03324           MSG("TrackSR",Msg::kDebug) <<  " track 2 " << track2->GetCluster(0)->GetPlane() << "/" << 
03325                     track2->GetCluster(track2->GetLast())->GetPlane() << endl;
03326             if(track2->GetLast()>fParmTrk2DNSeed && !track2->IsBad()){
03327               if(track1->GetPlaneView()==track2->GetPlaneView()){
03328                 merge=false;
03329                 if(track1!=track2){
03330                   if(direction>0){
03331                     fwd1=min(nremove1,track1->GetLast()-2);
03332                     bck1=max(track1->GetLast()-nremove1,2);
03333                     fwd2=min(nremove2,track2->GetLast()-2);
03334                     bck2=max(track2->GetLast()-nremove2,2);
03335                   }
03336                   else{
03337                     fwd1=max(track1->GetLast()-nremove1,2);
03338                     bck1=min(nremove1,track1->GetLast()-2);
03339                     fwd2=max(track2->GetLast()-nremove2,2);
03340                     bck2=min(nremove2,track2->GetLast()-2);
03341                   }
03342                   MSG("TrackSR",Msg::kDebug) << " Tracks " << itrack1 << "/" << itrack2 << " track 1 " << track1->GetCluster(fwd1)->GetPlane() << "/" << 
03343                     track1->GetCluster(bck1)->GetPlane() << "    track 2 " << track2->GetCluster(fwd2)->GetPlane() << "/" << track2->GetCluster(bck2)->GetPlane() << endl;
03344                   if(track2->GetCluster(bck2)->GetPlane()<
03345                      track1->GetCluster(fwd1)->GetPlane()){  // track1 downstream of track2
03346                     MSG("TrackSR",Msg::kDebug) << "track 1 is downstream of track 2" << endl;
03347                     
03348                     
03349                     Float_t upos,vpos,uslope,vslope;
03350                     if(track1->GetPlaneView()==PlaneView::kU){
03351                       upos = track1->GetCluster(fwd1)->GetTPos();
03352                       vpos = fHoughVIntercept + track1->GetCluster(fwd1)->GetZPos()*fHoughVSlope;
03353                       uslope=track1->GetForwardSlope(fwd1);
03354                       vslope=fHoughVSlope;
03355                     }
03356                     else{
03357                       vpos = track1->GetCluster(fwd1)->GetTPos();
03358                       upos = fHoughUIntercept + track1->GetCluster(fwd1)->GetZPos()*fHoughUSlope;
03359                       vslope=track1->GetForwardSlope(fwd1);
03360                       uslope=fHoughUSlope;
03361                     }
03362                     Int_t dir=-1;       
03363                     Int_t numSkipped=FindNumSkippedPlanes(track2->GetCluster(bck2)->GetPlane(),
03364                                                           track1->GetCluster(fwd1)->GetPlane(),
03365                                                           upos,vpos,
03366                                                           track1->GetCluster(fwd1)->GetZPos(),
03367                                                           uslope,
03368                                                           vslope,
03369                                                           dir,
03370                                                           track1->GetPlaneView());
03371                     MSG("TrackSR",Msg::kDebug) << "# skipped planes = " << numSkipped << " tpos " << track1->GetCluster(fwd1)->GetTPos() << " " << track2->GetCluster(bck2)->GetTPos() << endl;
03372                     if((TMath::Abs(track1->GetCluster(fwd1)->GetTPos())<0.375 && TMath::Abs(track2->GetCluster(bck2)->GetTPos())<0.375  ) || numSkipped<=fParmTrk2DNSkipHit ){ 
03373 
03374                       Float_t midZ=0.5*(track1->GetCluster(fwd1)->GetZPos()+
03375                                         track2->GetCluster(bck2)->GetZPos());
03376                       Float_t delZ=TMath::Abs(track1->GetCluster(fwd1)->GetZPos()-track2->GetCluster(bck2)->GetZPos());
03377                       Float_t textrap1=track1->GetCluster(fwd1)->GetTPos()+
03378                         track1->GetForwardSlope(fwd1)*(midZ-track1->GetCluster(fwd1)->GetZPos());        
03379                       Float_t textrap2=track2->GetCluster(bck2)->GetTPos()+
03380                         track2->GetForwardSlope(bck2)*(midZ-track2->GetCluster(bck2)->GetZPos());
03381                       Float_t delTPos=TMath::Abs(textrap1-textrap2);
03382                       Float_t tPosErr=  TMath::Sqrt(track1->GetCluster(fwd1)->GetTPosError()*track1->GetCluster(fwd1)->GetTPosError()+
03383                                                     track2->GetCluster(bck2)->GetTPosError()*track2->GetCluster(bck2)->GetTPosError());   
03384                       Float_t delSlope= TMath::Abs(track1->GetForwardSlope(fwd1)-track2->GetForwardSlope(bck2));
03385                       // multiply maxResid, maxSlopeDiff by root(2) since we are extrapolating two tracks to center, and not accounting for bending, mcs, etc 
03386                       Float_t maxResid = 1.4*(fParmTrk2DLinA+fParmTrk2DLinB*0.5*delZ);
03387                       Float_t maxSlopeDiff = 1.4*(TMath::Sqrt(tPosErr*tPosErr/(delZ*delZ)+fParmTrk2DDirCosError2));
03388                       MSG("TrackSR",Msg::kDebug) << " delTPos " << delTPos << "/" << maxResid << " delSlope " << delSlope <<"/" << maxSlopeDiff << endl;
03389                       if(delTPos<maxResid   &&
03390                          delSlope< maxSlopeDiff ){
03391                         merge=true;
03392                         MSG("TrackSR",Msg::kDebug) << " merging tracks " << endl;
03393                       }
03394                       
03395                     }
03396                     if(merge){
03397                       if( nremove1 || nremove2 )  MSG("TrackSR",Msg::kDebug) << " removing hits " << nremove1 << " " << nremove2 <<  endl;
03398                       for(Int_t i=0;i<nremove1;i++){
03399                         if(direction>0){
03400                           track1->RemoveAt(0);
03401                         }
03402                         else{
03403                           track1->RemoveAt(track1->GetLast());
03404                         }
03405                         track1->Compress();
03406                       }
03407                       for(Int_t i=0;i<nremove2;i++){
03408                         if(direction>0){
03409                           track2->RemoveAt(track2->GetLast());
03410                         }
03411                         else{
03412                           track2->RemoveAt(0);
03413                         }
03414                         track2->Compress();
03415                       }
03416                     }
03417                   }
03418                   else if (track1->GetCluster(bck1)->GetPlane()<
03419                            track2->GetCluster(fwd2)->GetPlane()){ // track2 downstream of track 1
03420                     MSG("TrackSR",Msg::kDebug) << "track 2 is downstream of track 1" << endl;
03421                     
03422                     Int_t dir=-1;
03423                     Float_t upos,vpos,uslope,vslope;
03424                     if(track1->GetPlaneView()==PlaneView::kU){
03425                       upos = track2->GetCluster(fwd2)->GetTPos();
03426                       vpos = fHoughVIntercept + track2->GetCluster(fwd2)->GetZPos()*fHoughVSlope;
03427                       uslope=track2->GetForwardSlope(fwd2);
03428                       vslope=fHoughVSlope;
03429                     }
03430                     else{
03431                       vpos = track2->GetCluster(fwd2)->GetTPos();
03432                       upos = fHoughUIntercept + track2->GetCluster(fwd2)->GetZPos()*fHoughUSlope;
03433                       vslope=track2->GetForwardSlope(fwd2);
03434                       uslope=fHoughUSlope;
03435                     }
03436                     Int_t numSkipped=FindNumSkippedPlanes(track1->GetCluster(bck1)->GetPlane(),
03437                                                           track2->GetCluster(fwd2)->GetPlane(),
03438                                                           upos,vpos,
03439                                                           track2->GetCluster(fwd2)->GetZPos(),
03440                                                           uslope,
03441                                                           vslope,
03442                                                           dir,
03443                                                           track1->GetPlaneView());
03444                     MSG("TrackSR",Msg::kDebug) << "# skipped planes = " << numSkipped << " tpos " << track1->GetCluster(bck1)->GetTPos() << " " << track2->GetCluster(fwd2)->GetTPos() << endl; 
03445                     if((TMath::Abs(track1->GetCluster(bck1)->GetTPos())<0.375 && TMath::Abs(track2->GetCluster(fwd2)->GetTPos())<0.375  ) || numSkipped<=fParmTrk2DNSkipHit){             
03446                       Float_t midZ=0.5*(track1->GetCluster(bck1)->GetZPos()+
03447                                         track2->GetCluster(fwd2)->GetZPos());
03448                       Float_t delZ=TMath::Abs(track1->GetCluster(bck1)->GetZPos()-track2->GetCluster(fwd2)->GetZPos());
03449                       Float_t textrap1=track1->GetCluster(bck1)->GetTPos()+
03450                         track1->GetForwardSlope(bck1)*(midZ-track1->GetCluster(bck1)->GetZPos());        
03451                       Float_t textrap2=track2->GetCluster(fwd2)->GetTPos()+
03452                         track2->GetForwardSlope(fwd2)*(midZ-track2->GetCluster(fwd2)->GetZPos());
03453                       Float_t delTPos = TMath::Abs(textrap1-textrap2);
03454                       Float_t tPosErr =  TMath::Sqrt(track1->GetCluster(bck1)->GetTPosError()*track1->GetCluster(bck1)->GetTPosError()+
03455                                                      track2->GetCluster(fwd2)->GetTPosError()*track2->GetCluster(fwd2)->GetTPosError());          
03456                       Float_t delSlope = TMath::Abs(track1->GetForwardSlope(bck1)-track2->GetForwardSlope(fwd2));
03457                       Float_t maxResid = 1.4*(fParmTrk2DLinA+fParmTrk2DLinB*0.5*delZ);
03458                       Float_t maxSlopeDiff = 1.4*(TMath::Sqrt(tPosErr*tPosErr/(delZ*delZ)+fParmTrk2DDirCosError2));
03459                       MSG("TrackSR",Msg::kDebug) << " delTPos " << delTPos << "/" << maxResid << " delSlope " << delSlope <<"/" << maxSlopeDiff << endl;
03460                       if(delTPos<maxResid   &&
03461                          delSlope< maxSlopeDiff && 
03462                          track1->GetLast()-nremove1>3 && track2->GetLast()-nremove2>3){
03463                         merge=true;
03464                         MSG("TrackSR",Msg::kDebug) << " merging tracks " << endl;
03465                       }         
03466                     }
03467                     if(merge){
03468                       if( nremove1 || nremove2 )  MSG("AlgTrackSRList",Msg::kDebug) << " removing hits " << nremove1 << " " << nremove2 <<  endl;
03469                       for(Int_t i=0;i<nremove1;i++){
03470                         if(direction>0){
03471                           track1->RemoveAt(track1->GetLast());
03472                         }
03473                         else{
03474                           track1->RemoveAt(0);
03475                         }
03476                         track1->Compress();
03477                       }
03478                       for(Int_t i=0;i<nremove2;i++){
03479                         if(direction>0){
03480                           track2->RemoveAt(0);
03481                         }
03482                         else{
03483                           track2->RemoveAt(track2->GetLast());
03484                         }
03485                         track2->Compress();
03486                       }
03487                     }     
03488                   }
03489                   if(merge){
03490                     for (Int_t icluster=0;icluster<=track2->GetLast();icluster++){
03491                       TrackClusterSR * tc0=track2->GetCluster(icluster);
03492                       Double_t slope = track2->GetForwardSlope(icluster); 
03493                       track1->Add(tc0,slope);
03494                     }
03495                     track2->SetIsBad(true);
03496                   }
03497                 }
03498               }
03499             }
03500             itrack2++;
03501           }
03502         }
03503         itrack1++;
03504       }
03505     }
03506   }
03507     //remove the bad tracks from the list
03508   for (Int_t itrk=0; itrk<=tracklist->GetLast(); itrk++) {
03509     track1 = dynamic_cast<Track2DSR*>(tracklist->At(itrk));
03510     if(track1->IsBad()){
03511       MSG("TrackSR" , Msg::kDebug) << "removing bad track " << itrk 
03512                                    << " with beg plane = "
03513                                    << track1->GetBegPlane() << endl;
03514       
03515       tracklist->RemoveAt(itrk);
03516       delete track1;
03517     }
03518   }
03519   tracklist->Compress();
03520 }

Bool_t AlgTrackSRList::PlaneIsActive ( Int_t  plane,
Float_t  u,
Float_t  v,
Float_t  tol 
) [private]

Definition at line 2413 of file AlgTrackSRList.cxx.

References PlaneOutline::DistanceToEdge(), fDetector, fPL, fVldc, PlexPlaneId::GetPlaneCoverage(), PlexPlaneId::GetPlaneView(), PlaneOutline::IsInside(), PlexPlaneId::IsValid(), and PlaneCoverage::kNoActive.

Referenced by FindNumSkippedPlanes().

02414 {
02415   UgliGeomHandle ugh(fVldc); 
02416   UgliScintPlnHandle planeid;  
02417   PlexPlaneId plexPlane(fDetector,plane, 0); 
02418   if(!plexPlane.IsValid()) return false; 
02419   if(plexPlane.GetPlaneCoverage() == PlaneCoverage::kNoActive) return false;
02420   if(projErr<0.3)projErr=0.3;
02421   float x = 0.707*(u-v);
02422   float y = 0.707*(u+v);
02423   float dist,xedge,yedge;
02424   fPL.DistanceToEdge(x, y,
02425                     plexPlane.GetPlaneView(),
02426                     plexPlane.GetPlaneCoverage(),
02427                     dist, xedge, yedge);
02428   Bool_t isInside = fPL.IsInside( x, y,
02429                     plexPlane.GetPlaneView(),
02430                    plexPlane.GetPlaneCoverage());
02431   isInside &= dist>projErr;
02432   return isInside;
02433        
02434 }

void AlgTrackSRList::RemoveStripsInSlice ( CandTrackSRListHandle cth,
CandSliceListHandle slicelist 
) [private]

Definition at line 1115 of file AlgTrackSRList.cxx.

References CandHandle::AddDaughterLink(), fDetector, CandHandle::GetDaughterIterator(), TrackClusterSR::GetPlane(), TrackClusterSR::GetStripList(), CandTrackSRListHandle::GetTrackClusterList(), CandHandle::GetUidInt(), CandHandle::IsEqual(), Detector::kNear, CandHandle::RemoveDaughter(), CandRecoHandle::SetCandSlice(), and tc.

Referenced by RunAlg().

01117 {
01118  
01119   if(!fDetector==Detector::kNear)return;
01120   if(slicelist){
01121     CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
01122     for (CandSliceHandle *slice = sliceItr(); slice ; slice = sliceItr()) {
01123       CandSliceHandle *  dupSlice=slice->DupHandle(); 
01124       CandStripHandleItr slicestripItr(dupSlice->GetDaughterIterator());
01125       for (CandStripHandle *strip = slicestripItr(); strip ; strip = slicestripItr()) {
01126         if(strip->GetPlane()>120){
01127           CandDigitHandleItr digitItr1(strip->GetDaughterIterator());
01128           CandDigitHandle *digit1 = dynamic_cast<CandDigitHandle*>(digitItr1());
01129           Bool_t stripIsOnTrack=false;
01130           Bool_t samePixel=false;  
01131 
01132         TObjArray *tclist = cth.GetTrackClusterList();
01133         for (int i=0; i<=tclist->GetLast(); i++) {
01134           TrackClusterSR *tc = dynamic_cast<TrackClusterSR*>(tclist->At(i));
01135           if(tc->GetPlane()>120 ){
01136             const TObjArray *trackstriplist = tc->GetStripList();
01137             for (int j=0; j<=trackstriplist->GetLast(); j++) {
01138               CandStripHandle * trackstrip = dynamic_cast<CandStripHandle*>(trackstriplist->At(j));
01139 
01140               CandDigitHandleItr digitItr2(trackstrip->GetDaughterIterator());
01141               CandDigitHandle *digit2 = dynamic_cast<CandDigitHandle*>(digitItr2());
01142 
01143               if(trackstrip->IsEqual(strip)){
01144                 stripIsOnTrack=true;
01145               }
01146               if(digit1->GetChannelId()==digit2->GetChannelId()){
01147                 samePixel=true;
01148               }
01149             }
01150           }
01151         }
01152         if(!stripIsOnTrack && samePixel){
01153           dupSlice->RemoveDaughter(strip);
01154           
01155         }
01156         }
01157       }
01158       if(dupSlice->GetUidInt()!=slice->GetUidInt()){
01159         CandTrackSRHandleItr trackItr(cth.GetDaughterIterator());
01160         for (CandTrackSRHandle *track = trackItr(); track ; track = trackItr()) {
01161           if(track->GetCandSlice()->GetUidInt()==slice->GetUidInt()){
01162             CandTrackSRHandle * newtrack = track->DupHandle();
01163             newtrack->SetCandSlice(dupSlice);
01164             cth.AddDaughterLink(*newtrack);
01165             cth.RemoveDaughter(track);
01166             delete newtrack;
01167           }
01168         }
01169         slicelist->AddDaughterLink(*dupSlice);
01170         slicelist->RemoveDaughter(slice);
01171       }
01172       delete dupSlice;
01173     } 
01174   }
01175 }

Int_t AlgTrackSRList::RemoveTrackSubsets ( Track2DSRItr &  trackItr1,
Track2DSRItr &  trackItr2,
TObjArray *  trackList 
) [private]

Definition at line 3138 of file AlgTrackSRList.cxx.

References Track2DSR::Compress(), Track2DSR::fIterate, fParmTrk2DSubsetNHit, fParmTrk2DTwinMatchFraction, Track2DSR::GetBegPlane(), Track2DSR::GetChi2(), Track2DSR::GetCluster(), Track2DSR::GetEndPlane(), Track2DSR::GetLast(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), Track2DSR::IsBad(), Msg::kDebug, MSG, and Track2DSR::SetIsBad().

Referenced by RunAlg().

03141 {
03142   Track2DSR *track1 = 0;
03143   Track2DSR *track2 = 0;
03144 
03145   trackItr1.ResetFirst();
03146   trackItr2.ResetFirst();
03147 
03148   Int_t removeCtr = 0;
03149 
03150   //counters to keep track of which track we are on in 
03151   //each set - dont want to remove a track because we happen to
03152   //be on the same one in the different sets
03153 
03154   Int_t itrk1 = 0;
03155   Int_t itrk2 = 0;
03156  
03157   Double_t chi21 = 0.;
03158   Double_t chi22 = 0.;
03159   
03160   Int_t nmatch = 0;
03161 
03162   TrackClusterSR *tc1 = 0;
03163   TrackClusterSR *tc2 = 0;
03164 
03165   Bool_t found = true;
03166   Bool_t match = false;
03167   Double_t matchfraction = 0.;
03168 
03169   //loop over tracks in both iterators to see if any match 
03170   while( (track1 = trackItr1()) ){
03171     
03172     //if there are any clusters in this track and it isnt a bad track
03173     if (track1->GetLast()>=0 && !track1->IsBad() ) {
03174 
03175        MSG("TrackSR", Msg::kDebug) << "track 1 " << "begPlane = " 
03176                                    << track1->GetBegPlane()
03177                                    << " end plane = " << track1->GetEndPlane() 
03178                                    << " itrk1  " << itrk1 << " last is  " 
03179                                    << track1->GetLast()
03180                                    << endl;
03181 
03182 
03183       itrk2 = 0;
03184       while( (track2 = trackItr2()) ){
03185     
03186         MSG("TrackSR", Msg::kDebug) << "track 2 " << "begPlane = " 
03187        
03188                                    << track2->GetBegPlane()
03189                                    << " end plane = " << track2->GetEndPlane() 
03190                                    << " itrk2 " << itrk2 << " last is  " << track2->GetLast() 
03191                                    << endl;
03192       
03193         if(track2->GetLast()>=0 && !track2->IsBad()){
03194           //if not on the same track but track2 has fewer clusters than
03195           //track1
03196           //if track 2 begins before track1 ends
03197           if( itrk1!=itrk2 
03198               && track2->GetLast()<=track1->GetLast() ){
03199                     
03200             MSG("TrackSR", Msg::kDebug) << "track2 could be a subset of track 1" 
03201                                         << endl; 
03202                 
03203             // check if track2 is subset of track1
03204             found = true;
03205             nmatch=0;
03206             
03207             //loop over the clusters in the 2 tracks to see if any are the same
03208             for(Int_t itc2=0; itc2<=track2->GetLast(); ++itc2){
03209               tc2 = track2->GetCluster(itc2);
03210               match = false;
03211               for(Int_t itc1=0; itc1<=track1->GetLast() && !match; ++itc1){
03212                 tc1 = track1->GetCluster(itc1);
03213                 if(tc1->GetPlane()==tc2->GetPlane() 
03214                    && tc1->GetMinStrip()==tc2->GetMinStrip() ){
03215                   match = true;
03216                   MSG("TrackSR" , Msg::kDebug) << "matching cluster found" << endl; 
03217                 }
03218               }//end loop over track1 clusters
03219 
03220               //if none are similar, no subset found
03221               if(!match){
03222                 MSG("TrackSR" , Msg::kDebug) << "matching cluster not found" 
03223                                              << endl;
03224                 found = false;
03225               }
03226               else ++nmatch;
03227               
03228             }//end loop over track2 clusters
03229                   
03230             //if found is still set then all the clusters in track2 were 
03231             //in track 1, so set track2 as bad
03232             if(found){
03233               track2->SetIsBad(true);
03234               MSG("TrackSR", Msg::kDebug) << "all track2 clusters found in track1 "
03235                                           << "set track " << itrk2 << " bad" 
03236                                           << endl;
03237             }
03238             else{
03239               matchfraction = (1.*nmatch)/(1.*track2->GetLast()+1.);
03240              MSG("TrackSR" , Msg::kDebug) << "matchfraction = " << matchfraction 
03241                                           << " twin match fraction = " 
03242                                           << fParmTrk2DTwinMatchFraction 
03243                                           << " nmatch = " << nmatch 
03244                                           << " subset n hits = " << fParmTrk2DSubsetNHit 
03245                                           << endl; 
03246 
03247               if(track2->GetLast()+1-nmatch<fParmTrk2DSubsetNHit && nmatch>0){
03248                 track2->SetIsBad(true);
03249                 MSG("TrackSR", Msg::kDebug) << "setting track2 bad" << endl;
03250                 
03251               }
03252               else if(matchfraction>fParmTrk2DTwinMatchFraction){
03253                 chi21 = track1->GetChi2()/(1.*track1->GetLast()+1.);
03254                 chi22 = track2->GetChi2()/(1.*track2->GetLast()+1.);
03255                 MSG("TrackSR", Msg::kDebug) << "chi21 = " << chi21
03256             
03257                                            << " chi22 = " << chi22 << endl;
03258                 if(chi21>chi22){
03259                   track1->SetIsBad(true);
03260                   MSG("TrackSR", Msg::kDebug) << "setting track1 bad, chi2" << endl;
03261                   
03262                 }
03263                 else{
03264                   track2->SetIsBad(true);
03265                   MSG("TrackSR", Msg::kDebug) << "setting track2 bad, chi2" << endl;
03266                  
03267                 }
03268               }//end if matchfraction makes it look like a twin
03269             }//end else for no matches found
03270           }//end if not the same track
03271         }//end if track2 is ok to use
03272         ++itrk2;
03273       }//end loop over trackItr2
03274       trackItr2.ResetFirst();
03275     }//end if track1 is ok to use
03276     ++itrk1;
03277   }//end loop over trackItr1
03278 
03279 
03280   //remove the bad tracks from the list
03281   for (Int_t itrk=0; itrk<=trackList->GetLast(); itrk++) {
03282     track1 = dynamic_cast<Track2DSR*>(trackList->At(itrk));
03283     if(track1->IsBad() || track1->fIterate<0){
03284 
03285       MSG("TrackSR" , Msg::kDebug) << "removing bad track " << itrk 
03286   
03287                                   << " with beg plane = "
03288                                   << track1->GetBegPlane() 
03289                                   << " is bad " << (Int_t)track1->IsBad()
03290                                   << " iteration " << track1->fIterate << endl;
03291 
03292       trackList->RemoveAt(itrk);
03293       delete track1;
03294       ++removeCtr;
03295     }
03296   }
03297   trackList->Compress();
03298 
03299   return removeCtr;
03300 }

void AlgTrackSRList::RemoveUnusedSpectStrips ( CandTrackSRListHandle cth,
CandStripListHandle striplist 
) [private]

Definition at line 1074 of file AlgTrackSRList.cxx.

References fDetector, CandHandle::GetDaughterIterator(), TrackClusterSR::GetPlane(), TrackClusterSR::GetStripList(), CandTrackSRListHandle::GetTrackClusterList(), CandHandle::IsEqual(), Detector::kNear, CandHandle::RemoveDaughter(), and tc.

Referenced by RunAlg().

01076 {
01077   if(!fDetector==Detector::kNear)return;
01078   if(striplist){
01079     CandStripHandleItr stripItr(striplist->GetDaughterIterator());
01080     for (CandStripHandle *strip = stripItr(); strip ; strip = stripItr()) {
01081       if(strip->GetPlane()>120){
01082         CandDigitHandleItr digitItr1(strip->GetDaughterIterator());
01083         CandDigitHandle *digit1 = dynamic_cast<CandDigitHandle*>(digitItr1());
01084         Bool_t stripIsOnTrack=false;
01085         Bool_t samePixel=false;
01086         TObjArray *tclist = cth.GetTrackClusterList();
01087         for (int i=0; i<=tclist->GetLast(); i++) {
01088           TrackClusterSR *tc = dynamic_cast<TrackClusterSR*>(tclist->At(i));
01089           if(tc->GetPlane()>120 ){
01090             const TObjArray *trackstriplist = tc->GetStripList();
01091             for (int j=0; j<=trackstriplist->GetLast(); j++) {
01092               CandStripHandle * trackstrip = dynamic_cast<CandStripHandle*>(trackstriplist->At(j));
01093               CandDigitHandleItr digitItr2(trackstrip->GetDaughterIterator());
01094               CandDigitHandle *digit2 = dynamic_cast<CandDigitHandle*>(digitItr2());
01095 
01096               if(trackstrip->IsEqual(strip)){
01097                 stripIsOnTrack=true;
01098               }
01099               if(digit1->GetChannelId()==digit2->GetChannelId()){
01100                 samePixel=true;
01101               }
01102             }
01103           }
01104         }
01105         if(!stripIsOnTrack && samePixel){
01106           striplist->RemoveDaughter(strip);
01107           
01108         }
01109       }
01110     }
01111   }
01112 }

void AlgTrackSRList::RunAlg ( AlgConfig ac,
CandHandle ch,
CandContext cx 
) [virtual]

Implements AlgBase.

Definition at line 189 of file AlgTrackSRList.cxx.

References AddClustersToTracks(), CleanUp(), DeleteTwinTracks(), fCDLHnew, fCSlcLHnew, fCSLHnew, fDetector, fHoughUIntercept, fHoughUSlope, fHoughVIntercept, fHoughVSlope, FillInGaps(), FillTrackList(), Find2DTrackEndPoints(), CandHandle::FindDaughter(), FindTimingDirection(), fMapIsWide, FormCandTrackSR(), fParmDiffViewBegPlaneMatch, fParmDiffViewEndPlaneMatch, fParmDiffViewPlaneOverlap, fParmDiffViewTimeMatch, fParmExtrapError, fParmHitNTime, fParmHitTime, fParmHoughMinBin, fParmHoughMinFrac, fParmHoughMinFracZoom, fParmHoughMinInterBinSize, fParmIsCosmic, fParmMaxNExtraStrip, fParmMaxNStrip, fParmMaxTimingResid, fParmMinClusterPulseHeight, fParmMinSingleHit, fParmMinStripPulseHeight, fParmMisalignmentError, fParmSingleHitDef, fParmTrk2DAlpha, fParmTrk2DDirCosError2, fParmTrk2DDirCosNSigma, fParmTrk2DEndPlaneDiff, fParmTrk2DHitFraction, fParmTrk2DLinA, fParmTrk2DLinA0, fParmTrk2DLinB, fParmTrk2DLinB0, fParmTrk2DMaxResid, fParmTrk2DNContiguous, fParmTrk2DNHough, fParmTrk2DNHough0, fParmTrk2DNSeed, fParmTrk2DNSigmaA, fParmTrk2DNSkip, fParmTrk2DNSkipHit, fParmTrk2DNSkipRemove, fParmTrk2DPlnEnd, fParmTrk2DSubsetDHit, fParmTrk2DSubsetNHit, fParmTrk2DTwinMatchFraction, fParmTrk2DWin0, fParmTrk3DTrackMax, fParmTrk3DTwinFrac, fParmTrkClsNSkip, fParmTrkNPlaneGoodDir, fParmUseShowerLikePlanes, fParmUseWideClusters, fSimFlag, fSliceDzDs, fSliceEndPlane, fSliceEndUPlane, fSliceEndVPlane, fSliceNPlanes, fSliceNUPlanes, fSliceNVPlanes, fSliceVtxPlane, fSliceVtxUPlane, fSliceVtxVPlane, fVldc, Registry::Get(), TrackClusterSR::GetBegTime(), CandRecord::GetCandHeader(), CandHandle::GetCandRecord(), CandContext::GetCandRecord(), Track2DSR::GetCluster(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), HoughViewSR::GetHoughTrack(), AlgFactory::GetInstance(), Registry::GetInt(), HoughTrackSR::GetLargestPeakIndex(), HoughViewSR::GetLargestTrackIndex(), Track2DSR::GetLast(), TrackClusterSR::GetMaxTPos(), TrackClusterSR::GetMinTPos(), CandContext::GetMom(), HoughViewSR::GetNTrack(), HoughTrackSR::GetPeakIntercept(), HoughTrackSR::GetPeakSlope(), TrackClusterSR::GetPlane(), Track2DSR::GetPlaneView(), VldContext::GetSimFlag(), CandHeader::GetSnarl(), TrackClusterSR::GetTPos(), CandHandle::GetUidInt(), RecMinos::GetVldContext(), CandHandle::GetVldContext(), TrackClusterSR::GetZPos(), IdentifyBadTracks(), Track2DSR::IsBad(), TrackClusterSR::IsValid(), TrackClusterSR::IsWide(), HoughViewSR::Iterate(), Msg::kDebug, Msg::kError, KeyOnClusterPlane(), KeyOnUCluster(), KeyOnUClusterNotSL(), KeyOnUClusterNotWide(), KeyOnUClusterNotWideSL(), KeyOnUViewTracks(), KeyOnUViewVetoCharge(), KeyOnUViewVetoChargeNear(), KeyOnVCluster(), KeyOnVClusterNotSL(), KeyOnVClusterNotWide(), KeyOnVClusterNotWideSL(), KeyOnVViewTracks(), KeyOnVViewVetoCharge(), KeyOnVViewVetoChargeNear(), Detector::kNear, PlaneView::kU, SimFlag::kUnknown, Detector::kUnknown, PlaneView::kV, MakeTrackClusters(), MarkUsedClusters(), MergeTracks(), Munits::mm, MSG, Munits::ns, RemoveStripsInSlice(), RemoveTrackSubsets(), RemoveUnusedSpectStrips(), CandRecord::SecureCandHandle(), HoughViewSR::SetClusterList(), HoughViewSR::SetMinInterBinSize(), CandHandle::SetName(), HoughViewSR::SetPeakMin(), HoughViewSR::SetPeakMinFrac(), HoughViewSR::SetPeakMinFracZoom(), SpectrometerTracking(), CandStripHandle::StripSRKeyFromBegTime(), CandStripHandle::StripSRKeyFromPSEId(), tc, and LinearFit::Weighted().

00190 {
00191   MSG("Alg", Msg::kDebug) << "Starting AlgTrackSRList::RunAlg()" << endl;
00192 
00193   assert(cx.GetDataIn());
00194 
00195   if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
00196     return;
00197   }
00198 
00199   CandTrackSRListHandle &cth = dynamic_cast<CandTrackSRListHandle &>(ch);
00200 
00201   const CandSliceListHandle *slicelist = 0;
00202   const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
00203   for (Int_t i=0; i<=cxin->GetLast(); i++) {
00204     TObject *tobj = cxin->At(i);
00205     if (tobj->InheritsFrom("CandSliceListHandle")) {
00206       slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
00207     }
00208   }
00209   fDetector = Detector::kUnknown;
00210   fSimFlag = SimFlag::kUnknown;
00211 
00212   if (!slicelist) {
00213     MSG("error",Msg::kError) <<
00214       "CandSliceListHandle missing\n";
00215   }
00216   else {
00217     fDetector = slicelist->GetVldContext()->GetDetector();
00218     fSimFlag = slicelist->GetVldContext()->GetSimFlag();
00219     MSG("TrackSR", Msg::kDebug) << "event number " 
00220                                  << slicelist->GetCandRecord()->GetCandHeader()->GetSnarl() 
00221                                  << endl;
00222   }
00223 
00224 
00225 // Create Candcontext
00226   CandContext cxx(this,cx.GetMom());
00227 
00228 // Get singleton instance of AlgFactory
00229   AlgFactory &af = AlgFactory::GetInstance();
00230 
00231   fParmMaxNStrip = ac.GetInt("MaxNStrip");
00232   fParmMaxNExtraStrip = ac.GetInt("MaxNExtraStrip");
00233   fParmTrk2DPlnEnd = ac.GetInt("Trk2DPlnEnd");
00234   fParmTrk2DWin0 = ac.GetDouble("Trk2DWin0");
00235   fParmHitNTime = ac.GetDouble("HitNTime");
00236   fParmHitTime = ac.GetDouble("HitTime");
00237   fParmTrk2DNSkip = ac.GetInt("Trk2DNSkip");
00238   fParmTrk2DNSkipHit = ac.GetInt("Trk2DNSkipHit");
00239   fParmTrk2DAlpha = ac.GetDouble("Trk2DAlpha");
00240   fParmTrk2DNSkipRemove = ac.GetInt("Trk2DNSkipRemove");
00241   fParmTrkNPlaneGoodDir = ac.GetInt("TrkNPlaneGoodDir");
00242   fParmTrk2DDirCosNSigma = ac.GetDouble("Trk2DDirCosNSig") ;
00243   fParmTrk2DDirCosError2 = ac.GetDouble("Trk2DDirCosError")*ac.GetDouble("Trk2DDirCosError");
00244   fParmTrk2DEndPlaneDiff = ac.GetInt("Trk2DEndPlaneDiff");
00245   fParmMinSingleHit = ac.GetInt("MinSingleHit");
00246   fParmSingleHitDef = ac.GetInt("SingleHitDef");
00247   fParmTrk2DNHough0 = ac.GetInt("Trk2DNHough0");
00248   fParmTrk2DLinA0 = ac.GetDouble("Trk2DLinA0");
00249   fParmTrk2DLinB0 = ac.GetDouble("Trk2DLinB0");
00250   fParmTrk2DNHough = ac.GetInt("Trk2DNHough");
00251   fParmTrk2DLinA = ac.GetDouble("Trk2DLinA");
00252   fParmTrk2DLinB = ac.GetDouble("Trk2DLinB");
00253   fParmTrk2DNSigmaA = ac.GetDouble("Trk2DNSigmaA");
00254   fParmTrk2DNSeed = ac.GetInt("Trk2DNSeed");
00255   fParmTrk2DMaxResid = ac.GetDouble("Trk2DMaxResid");
00256   fParmTrk2DNContiguous = ac.GetInt("Trk2DNContiguous");
00257   fParmTrk2DHitFraction = ac.GetDouble("Trk2DHitFraction");
00258   fParmTrk2DTwinMatchFraction = ac.GetDouble("Trk2DTwinMatchFraction");
00259   fParmTrk2DSubsetNHit = ac.GetInt("Trk2DSubsetNHit");
00260   fParmTrk2DSubsetDHit = ac.GetInt("Trk2DSubsetDHit");
00261   fParmDiffViewBegPlaneMatch = ac.GetInt("DiffViewBegPlaneMatch");
00262   fParmDiffViewEndPlaneMatch = ac.GetInt("DiffViewEndPlaneMatch");
00263   fParmDiffViewPlaneOverlap = ac.GetDouble("DiffViewPlaneOverlap");
00264   fParmDiffViewTimeMatch = ac.GetDouble("DiffViewTimeMatch");
00265   fParmTrk3DTwinFrac = ac.GetDouble("Trk3DTwinFrac");
00266   fParmTrkClsNSkip = ac.GetInt("TrkClsNSkip");
00267   fParmTrk3DTrackMax = ac.GetInt("Trk3DTrackMax");
00268   fParmIsCosmic = ac.GetInt("IsCosmic");
00269   fParmMinStripPulseHeight = ac.GetDouble("MinStripPulseHeight");
00270   fParmMinClusterPulseHeight = ac.GetDouble("MinClusterPulseHeight");
00271   fParmUseWideClusters = ac.GetInt("UseWideClusters");
00272   fParmUseShowerLikePlanes = ac.GetInt("UseShowerLikePlanes");
00273   fParmMaxTimingResid = ac.GetInt("MaxTimingResid")*Munits::ns;
00274   fParmMisalignmentError = ac.GetInt("MisalignmentError")*Munits::mm;
00275   fParmExtrapError = ac.GetDouble("ExtrapError");
00276 
00277   fParmHoughMinBin =  fParmTrk2DNHough0 ;
00278   fParmHoughMinFrac = 0.5;
00279   fParmHoughMinFracZoom = .75;
00280   fParmHoughMinInterBinSize = .02;
00281 
00282   int EnableSpectFind = 1;
00283   EnableSpectFind = ac.GetInt("EnableSpectFind");
00284 
00285   const CandRecord *candrec = cx.GetCandRecord();
00286   assert(candrec);
00287   const VldContext *vldcptr = candrec->GetVldContext();
00288   assert(vldcptr);
00289   fVldc = *vldcptr; 
00290   
00291   UgliGeomHandle ugh(fVldc);
00292 
00293  // Get an AlgHandle to AlgTrackSR with proper AlgConfig
00294   const char *pTrackAlgConfig = 0;
00295   ac.Get("TrackAlgConfig",pTrackAlgConfig);
00296   AlgHandle ah = af.GetAlgHandle("AlgTrackSR",pTrackAlgConfig);
00297 
00298 
00299   // set new strip and digit list handles to null at this point
00300   // if we have near detector spectrometer tracking these will get set 
00301   // during process of the first slice
00302   fCDLHnew=0;
00303   fCSLHnew=0;
00304   fCSlcLHnew=0;
00305 
00306   //get an iterator over the slice list passed into the algorithm
00307   //and loop over all slices
00308   Int_t islice=0;
00309 
00310   CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
00311   while (CandSliceHandle *slice = sliceItr()) {
00312     MSG("TrackSR", Msg::kDebug) << " ******* Slice ********* " << islice <<  endl;
00313   
00314     islice++;
00315 
00316     //make variables to keep track of the event extent in 
00317     //the different views.  also keep track of the first and 
00318     //last planes hit in the event and the number of planes 
00319     //hit in a given view.
00320     fSliceVtxUPlane = 1000;
00321     fSliceVtxVPlane = 1000;
00322     fSliceEndUPlane = -1;
00323     fSliceEndVPlane = -1;
00324     fSliceVtxPlane = 1000;
00325     fSliceEndPlane = -1;
00326     fSliceNUPlanes = 0;
00327     fSliceNVPlanes = 0;
00328     fSliceNPlanes = 0;
00329     
00330     //zero the map of wide planes for this slice
00331     for(Int_t i = 0; i<1000; ++i){
00332       fMapIsWide[i] = 0;
00333     }
00334     
00335     //get a couple of iterators over the strips in the slice    
00336     CandStripHandleItr uStripIter(slice->GetDaughterIterator());        
00337     CandStripHandleItr vStripIter(slice->GetDaughterIterator());        
00338     CandStripHandleItr allStripIter(slice->GetDaughterIterator());      
00339         
00340      //sort the strips using the 2 function sort in navigation  
00341         uStripIter.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime);
00342         vStripIter.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime);
00343         allStripIter.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime);
00344 
00345     /* for now leave in the old sort just in case       
00346     //set the key functions to sort the strips by plane, strip and time
00347     CandStripHandleKeyFunc *uStripKF = uStripIter.CreateKeyFunc();
00348     uStripKF->SetFun(CandStripHandle::KeyFromPlaneStripTime);
00349     uStripIter.GetSet()->AdoptSortKeyFunc(uStripKF);
00350     uStripKF = 0;
00351 
00352     CandStripHandleKeyFunc *vStripKF = vStripIter.CreateKeyFunc();
00353     vStripKF->SetFun(CandStripHandle::KeyFromPlaneStripTime); 
00354     vStripIter.GetSet()->AdoptSortKeyFunc(vStripKF);
00355     vStripKF = 0;
00356 
00357     CandStripHandleKeyFunc *allStripKF = allStripIter.CreateKeyFunc();
00358     allStripKF->SetFun(CandStripHandle::KeyFromPlaneStripTime);
00359     allStripIter.GetSet()->AdoptSortKeyFunc(allStripKF);
00360     allStripKF = 0;
00361     */
00362     allStripIter.ResetFirst();
00363    
00364     // a loop to check that the sorting worked
00365     /*
00366     CandStripHandle *strip = 0;
00367     while( (strip = allStripIter()) ){
00368       MSG("TrackSR", Msg::kDebug) << "plane " << strip->GetPlane() << " strip " 
00369                                   << strip->GetStrip() << " " 
00370                                   << (Int_t)strip->GetPlaneView() << " " 
00371                                   << (Int_t)PlaneView::kX << " " 
00372                                   << (Int_t)PlaneView::kY << " " 
00373                                   << (Int_t)PlaneView::kU << " " 
00374                                   << (Int_t)PlaneView::kV << " " 
00375                                   << endl;
00376     } 
00377     */
00378 
00379     //now adopt a select key function to only get those strips
00380     //in each view that were not vetoed 
00381     CandStripHandleKeyFunc *uSelectKF = uStripIter.CreateKeyFunc();
00382     if(fDetector!=Detector::kNear){
00383       uSelectKF->SetFun(KeyOnUViewVetoCharge);
00384     }
00385     else{
00386       uSelectKF->SetFun(KeyOnUViewVetoChargeNear);
00387     }
00388 
00389     uStripIter.GetSet()->AdoptSelectKeyFunc(uSelectKF);
00390     uSelectKF = 0;
00391 
00392     CandStripHandleKeyFunc *vSelectKF = vStripIter.CreateKeyFunc();
00393     if(fDetector!=Detector::kNear){
00394       vSelectKF->SetFun(KeyOnVViewVetoCharge);
00395     }
00396     else{
00397       vSelectKF->SetFun(KeyOnVViewVetoChargeNear);
00398     }
00399     vStripIter.GetSet()->AdoptSelectKeyFunc(vSelectKF);
00400     vSelectKF = 0;      
00401     
00402     uStripIter.ResetFirst();
00403     vStripIter.ResetFirst();
00404     allStripIter.ResetFirst();
00405 
00406     //find the vtx and end planes, as well
00407     //as the total number of planes in each view
00408     Find2DTrackEndPoints(uStripIter, 
00409                          fSliceVtxUPlane, 
00410                          fSliceEndUPlane, 
00411                          fSliceNUPlanes);
00412     Find2DTrackEndPoints(vStripIter, 
00413                          fSliceVtxVPlane, 
00414                          fSliceEndVPlane, 
00415                          fSliceNVPlanes);
00416     MSG("TrackSR", Msg::kDebug) << "planes in slice " << fSliceNVPlanes << " " 
00417                                 << fSliceNUPlanes << " " << fParmTrk2DNSeed << endl;
00418 
00419     if(fSliceNUPlanes<fParmTrk2DNSeed || fSliceNVPlanes<fParmTrk2DNSeed){
00420       MSG("TrackSR", Msg::kDebug) << "too few planes in slice" << endl;
00421       continue;
00422     }
00423 
00424     fSliceNPlanes = fSliceNUPlanes + fSliceNVPlanes;
00425 
00426     //get the rough tracking information from the HoughViewSR objects
00427     //these give you an idea of the slope and intercept for the 2D tracks
00428     //in the event, as well as the number of tracks to expect in the event
00429     //the SetStripList method returns true if there are enough planes in
00430     //each view that pass the requirements for estimating using the hough transform
00431 
00432     HoughViewSR houghUView;
00433     houghUView.SetPeakMin(fParmHoughMinBin);
00434     houghUView.SetPeakMinFrac(fParmHoughMinFrac);
00435     houghUView.SetPeakMinFracZoom(fParmHoughMinFracZoom);
00436     houghUView.SetMinInterBinSize(fParmHoughMinInterBinSize);
00437 
00438     HoughViewSR houghVView;
00439     houghVView.SetPeakMin(fParmHoughMinBin);
00440     houghVView.SetPeakMinFrac(fParmHoughMinFrac);
00441     houghVView.SetPeakMinFracZoom(fParmHoughMinFracZoom);
00442     houghVView.SetMinInterBinSize(fParmHoughMinInterBinSize);
00443 
00444     //quit trying to track if you dont have enough 
00445 
00446     Int_t maxUTrack = 0;
00447     Int_t maxVTrack = 0;
00448     
00449     Double_t uTrackSlope = 0.;
00450     Double_t uTrackIntercept = 0.;
00451     Double_t vTrackSlope = 0.;
00452     Double_t vTrackIntercept = 0.;
00453 
00454     //get the slopes and intercepts for the different views
00455  
00456     TObjArray *utrackClusterList = new TObjArray(1,0);   
00457     //fill the cluster list one view at a time
00458 
00459     MakeTrackClusters(utrackClusterList, uStripIter, cth, 1, 1);
00460     TrackClusterSR *tc = 0;
00461     Double_t xfit[500];
00462     Double_t yfit[500];
00463     Double_t wfit[500];
00464     Double_t parm[2];
00465     Double_t eparm[2] = {0.,0.};
00466     Double_t wtmp;
00467     Int_t ifit=0;
00468     for(Int_t i=0; i<=utrackClusterList->GetLast(); ++i){
00469       tc = dynamic_cast<TrackClusterSR*>(utrackClusterList->At(i));
00470       if(ifit<500){
00471         Bool_t use=true;
00472         for(Int_t j=0; j<=utrackClusterList->GetLast(); ++j){
00473           TrackClusterSR * tc2 = dynamic_cast<TrackClusterSR*>(utrackClusterList->At(j));
00474           if(i!=j  && tc->GetPlane()==tc2->GetPlane()){
00475             use=false;
00476             break;
00477           }
00478         }
00479         if(use){
00480           xfit[ifit] = tc->GetZPos();
00481           yfit[ifit] = tc->GetTPos();
00482           wtmp = 0.288*(tc->GetMaxTPos()-tc->GetMinTPos());
00483           wfit[ifit]=1.0/(wtmp*wtmp);
00484           ifit++;
00485         }
00486       }
00487     }
00488     if(ifit>2){
00489       LinearFit::Weighted(ifit,xfit,yfit,wfit,parm,eparm);
00490       uTrackSlope=parm[1];
00491       uTrackIntercept=parm[0];
00492     }
00493 
00494     houghUView.SetClusterList(utrackClusterList, fParmMinStripPulseHeight); 
00495 
00496     houghUView.Iterate();
00497     MSG("TrackSR", Msg::kDebug) << "u tracks " << houghUView.GetNTrack() << endl;
00498     if( houghUView.GetNTrack() > 0 ){
00499       maxUTrack = houghUView.GetNTrack();
00500       const HoughTrackSR *ht = houghUView.GetHoughTrack(houghUView.GetLargestTrackIndex());
00501       uTrackSlope = ht->GetPeakSlope(ht->GetLargestPeakIndex());
00502       uTrackIntercept = ht->GetPeakIntercept(ht->GetLargestPeakIndex());
00503     }
00504 
00505     TObjArray *vtrackClusterList = new TObjArray(1,0);   
00506     //fill the cluster list one view at a time
00507     MakeTrackClusters(vtrackClusterList, vStripIter, cth, 1, 1);
00508    
00509     ifit=0;
00510     for(Int_t i=0; i<=vtrackClusterList->GetLast(); ++i){
00511       tc = dynamic_cast<TrackClusterSR*>(vtrackClusterList->At(i));
00512       if(ifit<500){
00513         Bool_t use=true;
00514         for(Int_t j=0; j<=vtrackClusterList->GetLast(); ++j){
00515           TrackClusterSR * tc2 = dynamic_cast<TrackClusterSR*>(vtrackClusterList->At(j));
00516           if(i!=j  && tc->GetPlane()==tc2->GetPlane()){
00517             use=false;
00518             break;
00519           }
00520         }
00521         if(use){
00522           xfit[ifit] = tc->GetZPos();
00523           yfit[ifit] = tc->GetTPos();
00524           wtmp = 0.288*(tc->GetMaxTPos()-tc->GetMinTPos());
00525           wfit[ifit]=1.0/(wtmp*wtmp);
00526           ifit++;
00527         }
00528       }
00529     }
00530     if(ifit>2){
00531       LinearFit::Weighted(ifit,xfit,yfit,wfit,parm,eparm);
00532       vTrackSlope=parm[1];
00533       vTrackIntercept=parm[0];
00534     }
00535     houghVView.SetClusterList(vtrackClusterList, fParmMinStripPulseHeight); 
00536     houghVView.Iterate();
00537     MSG("TrackSR", Msg::kDebug) << "v tracks " << houghVView.GetNTrack() << endl;
00538     if( houghVView.GetNTrack() > 0 ){
00539       maxVTrack = houghVView.GetNTrack();
00540       const HoughTrackSR *ht = houghVView.GetHoughTrack(houghVView.GetLargestTrackIndex());
00541       vTrackSlope = ht->GetPeakSlope(ht->GetLargestPeakIndex());
00542       vTrackIntercept = ht->GetPeakIntercept(ht->GetLargestPeakIndex());
00543     }
00544 
00545 
00546  // delete track clusters
00547     for(Int_t i=0; i<=utrackClusterList->GetLast(); ++i){
00548       tc = dynamic_cast<TrackClusterSR*>(utrackClusterList->At(i));
00549       delete tc;
00550     }
00551     delete utrackClusterList;
00552 
00553     for(Int_t i=0; i<=vtrackClusterList->GetLast(); ++i){
00554       tc = dynamic_cast<TrackClusterSR*>(vtrackClusterList->At(i));
00555       delete tc;
00556     }
00557     delete vtrackClusterList;
00558     
00559     
00560     //get the slope of the track wrt the z direction, dz/ds
00561     fSliceDzDs = 1./TMath::Sqrt(1.+uTrackSlope*uTrackSlope+vTrackSlope*vTrackSlope);
00562 
00563     MSG("TrackSR",Msg::kDebug) << "Hough: # tracks = " << maxUTrack 
00564                                  << " " << maxVTrack 
00565                                  << endl
00566                                  << "Hough: slope = " << uTrackSlope 
00567                                  << " " << vTrackSlope 
00568                                  << endl
00569                                  << "Hough: intercept = " << uTrackIntercept << " " 
00570                                  << vTrackIntercept << endl;
00571     
00572     if(uTrackSlope==-1)uTrackSlope=0;
00573     if(vTrackSlope==-1)vTrackSlope=0;
00574 
00575     fHoughUSlope = uTrackSlope;
00576     fHoughVSlope = vTrackSlope;
00577     fHoughUIntercept = uTrackIntercept;
00578     fHoughVIntercept = vTrackIntercept;
00579 
00580 
00581     //set the initial fit direction with respect to z 
00582     //1 = forward, -1 = backwards
00583     Int_t idir = 1; 
00584     //if this is a beam file, just set the direction to -1 because
00585     //it is easier to work backwards to the vertex
00586     if(!fParmIsCosmic) idir = -1; 
00587 
00588     //figure out where the vtx and end planes for the event
00589     fSliceVtxPlane = TMath::Min(fSliceVtxUPlane,fSliceVtxVPlane);
00590     fSliceEndPlane = TMath::Max(fSliceEndUPlane,fSliceEndVPlane);
00591 
00592     //determine whether it is worth it to do the timing for the event
00593     //make sure at least 1 track from the hough transform in each view
00594     if(maxUTrack>0 && maxVTrack>0){
00595       allStripIter.ResetFirst();
00596       //check and see if you have the direction right using the
00597       //timing in the event
00598       if(FindTimingDirection(allStripIter, uTrackSlope, uTrackIntercept,
00599                              vTrackSlope, vTrackIntercept)<0 ){
00600         if(fParmIsCosmic){
00601           idir = -1;
00602           fSliceEndPlane = TMath::Min(fSliceVtxUPlane,fSliceVtxVPlane);
00603           fSliceVtxPlane = TMath::Max(fSliceEndUPlane,fSliceEndVPlane);
00604         }
00605       }
00606     }//end if enough tracks to do the timing
00607     
00608     allStripIter.ResetFirst();
00609 
00610     Int_t uExpectedStrips = fParmMaxNStrip-fParmMaxNExtraStrip;
00611     Int_t vExpectedStrips = fParmMaxNStrip-fParmMaxNExtraStrip;
00612 
00613     MSG("TrackSR", Msg::kDebug) << "max u expected strips " << uExpectedStrips
00614                                << " max v expected strips " << vExpectedStrips 
00615                                << endl;
00616 
00617     //the number of expected strips is based on the slope of the track in
00618     //tpos vs z /4.108 where 4.108 is the width of a strip in tpos
00619     //so the value is how many strips you think will be hit/plane
00620     Int_t expectedStrips = (Int_t)(TMath::Abs(uTrackSlope)/4.108)+1;
00621     
00622     //only change the # of expected strips if you did the hough tracking
00623     //and it returned with at least 1 track in each view
00624     if(expectedStrips<uExpectedStrips && maxUTrack>0) 
00625       uExpectedStrips = expectedStrips;
00626     
00627     expectedStrips = (Int_t)(TMath::Abs(vTrackSlope)/4.108)+1;
00628     
00629     if(expectedStrips<vExpectedStrips && maxVTrack>0) 
00630       vExpectedStrips = expectedStrips;
00631 
00632     MSG("TrackSR", Msg::kDebug) << "u expected strips " << uExpectedStrips
00633                                << " v expected strips " << vExpectedStrips 
00634                                << " u track slope " << uTrackSlope 
00635                                << " v track slope " << vTrackSlope << endl;
00636       
00637     //now we can make track clusters.  declare an array to hold the cluster
00638     //objects, but fill it in the MakeTrackClusters() method
00639     TObjArray *trackClusterList = new TObjArray(1,0);
00640    
00641     //fill the cluster list one view at a time
00642     MakeTrackClusters(trackClusterList, uStripIter, cth, uExpectedStrips, idir);
00643     MakeTrackClusters(trackClusterList, vStripIter, cth, vExpectedStrips, idir);
00644       
00645     //an iterator over all clusters in the event
00646     TrackClusterSRItr masterClusterItr(trackClusterList);
00647     TrackClusterSRKeyFunc *masterClusterKf = masterClusterItr.CreateKeyFunc();
00648     masterClusterKf->SetFun(KeyOnClusterPlane);
00649     masterClusterItr.GetSet()->AdoptSortKeyFunc(masterClusterKf);
00650     masterClusterKf = 0;
00651     
00652     //clusterItr has all clusters in a view that pass the given cuts
00653     TrackClusterSRItr clusterItr(trackClusterList);
00654     TrackClusterSRKeyFunc *clusterKf = clusterItr.CreateKeyFunc();
00655     clusterKf->SetFun(KeyOnClusterPlane);
00656     clusterItr.GetSet()->AdoptSortKeyFunc(clusterKf);
00657     clusterKf = 0;
00658  
00659  
00660     //planeClusterItr will keep track of the clusters in a given plane
00661     TrackClusterSRItr planeClusterItr(trackClusterList);
00662     TrackClusterSRKeyFunc *planeClusterKf = planeClusterItr.CreateKeyFunc();
00663     planeClusterKf->SetFun(KeyOnClusterPlane);
00664     planeClusterItr.GetSet()->AdoptSortKeyFunc(planeClusterKf);
00665     planeClusterKf = 0;
00666     
00667     //fill array of ints to ints where the index is the 
00668     //plane number and the value is 0 for planes
00669     //containing only good clusters, and 1 for planes with
00670     //wide clusters
00671     
00672     clusterItr.ResetFirst();
00673     TrackClusterSR *cluster = 0;
00674     while( (cluster = clusterItr()) ){
00675       if( cluster->IsWide() ) fMapIsWide[cluster->GetPlane()] = (Int_t)cluster->IsWide();
00676 
00677       MSG("TrackSR", Msg::kDebug) << "cluster on plane " << cluster->GetPlane() 
00678                                  << " tpos " << cluster->GetTPos() << " fMapIsWide "
00679                                  << fMapIsWide[cluster->GetPlane()] << endl;
00680     }
00681     clusterItr.ResetFirst();
00682     
00683     //check to see how many planes you might expect in a track
00684     //by looking for gaps in the planes hit.
00685  
00686     Int_t nmax3dtrk = 0;
00687     int trk2dmin = fParmTrk2DNSeed;
00688    
00689     
00691     //
00692     // 2D tracking
00693     //
00695     
00696     
00697     //loop over each view separately to do the 2d tracking
00698     Int_t viewCtr = 0;
00699     Double_t houghSlope = 0.;
00700     Double_t houghIntercept = 0.;
00701     Int_t nmaxtrack = 0;
00702     Int_t sliceNPlane = 0;
00703     Int_t sliceVtxPlane = 0;
00704     Int_t sliceEndPlane = 0;
00705     
00706     //create a TObjArray to hold the 2d tracks we find
00707     TObjArray *trackList = new TObjArray(1,0);
00708     TObjArray *seedClusters = new TObjArray(1,0);
00709     
00710     while(viewCtr<2){
00711       
00712       //first we want to select only those clusters for this view
00713       //we dont want to use any wide clusters when doing the track finding
00714       //depending on whether we will allow clusters from showerlike planes
00715       //we either use KeyOnXClusterNotWide or KeyOnXClusterNotWideSL
00716       //where the X is either U or V and NotWide excludes only wide
00717       //clusters, NotWideSL excludes wide clusters and those from showerlike
00718       //planes.  if we allow wide clusters, then use KeyOnXCluster to just
00719       //select a view
00720       
00721       //the master cluster list has all clusters that pass the min 
00722       //pulseheight cut.  if wide clusters are allowed it will have 
00723       //those too.  it will always have clusters from showerlike planes.
00724       //so it either uses KeyOnU(V)Cluster or KeyOnU(V)ClusterNotWide
00725       
00726       TrackClusterSRKeyFunc *selectKF = clusterItr.CreateKeyFunc();
00727 
00728       TrackClusterSRKeyFunc *planeSelectKF = planeClusterItr.CreateKeyFunc();
00729       TrackClusterSRKeyFunc *masterSelectKF = masterClusterItr.CreateKeyFunc();
00730 
00731       //u view first
00732       if(viewCtr == 0){
00733         houghSlope = uTrackSlope;
00734         houghIntercept = uTrackIntercept;
00735         nmaxtrack = maxUTrack;
00736         sliceNPlane = fSliceNUPlanes;
00737         sliceVtxPlane = fSliceVtxUPlane;
00738         sliceEndPlane = fSliceEndUPlane;
00739         
00740         if(fParmUseWideClusters) masterSelectKF->SetFun(KeyOnUCluster);
00741         else masterSelectKF->SetFun(KeyOnUClusterNotWide);
00742         
00743         if(fParmUseWideClusters && fParmUseShowerLikePlanes){
00744           selectKF->SetFun(KeyOnUCluster);
00745           planeSelectKF->SetFun(KeyOnUCluster);
00746         }
00747         else if(fParmUseWideClusters && !fParmUseShowerLikePlanes){
00748           selectKF->SetFun(KeyOnUClusterNotSL);
00749           planeSelectKF->SetFun(KeyOnUClusterNotSL);
00750         }
00751         else if(!fParmUseWideClusters && fParmUseShowerLikePlanes){
00752           selectKF->SetFun(KeyOnUClusterNotWide);
00753           planeSelectKF->SetFun(KeyOnUClusterNotWide);
00754         }
00755         else if(!fParmUseWideClusters && !fParmUseShowerLikePlanes){
00756           selectKF->SetFun(KeyOnUClusterNotWideSL);
00757           planeSelectKF->SetFun(KeyOnUClusterNotWideSL);
00758         }
00759       }//end if 0 view
00760       else if( viewCtr == 1){
00761         houghSlope = vTrackSlope;
00762         houghIntercept = vTrackIntercept;
00763         nmaxtrack = maxVTrack;
00764         sliceNPlane = fSliceNVPlanes;
00765         sliceVtxPlane = fSliceVtxVPlane;
00766         sliceEndPlane = fSliceEndVPlane;
00767         
00768         if(fParmUseWideClusters) masterSelectKF->SetFun(KeyOnVCluster);
00769         else masterSelectKF->SetFun(KeyOnVClusterNotWide);
00770         
00771         if(fParmUseWideClusters && fParmUseShowerLikePlanes){
00772           selectKF->SetFun(KeyOnVCluster);
00773           planeSelectKF->SetFun(KeyOnVCluster);
00774         }
00775         else if(fParmUseWideClusters && !fParmUseShowerLikePlanes){
00776           selectKF->SetFun(KeyOnVClusterNotSL);
00777           planeSelectKF->SetFun(KeyOnVClusterNotSL);
00778         }
00779         else if(!fParmUseWideClusters && fParmUseShowerLikePlanes){
00780           selectKF->SetFun(KeyOnVClusterNotWide);
00781           planeSelectKF->SetFun(KeyOnVClusterNotWide);
00782         }
00783         else if(!fParmUseWideClusters && !fParmUseShowerLikePlanes){
00784           selectKF->SetFun(KeyOnVClusterNotWideSL);
00785             planeSelectKF->SetFun(KeyOnVClusterNotWideSL);
00786         }
00787       }//end if v view
00788       
00789       clusterItr.GetSet()->AdoptSelectKeyFunc(selectKF);
00790       planeClusterItr.GetSet()->AdoptSelectKeyFunc(planeSelectKF);
00791       masterClusterItr.GetSet()->AdoptSelectKeyFunc(masterSelectKF);
00792       selectKF = 0;
00793       planeSelectKF = 0;
00794       masterSelectKF = 0;
00795       
00796       TrackClusterSR *cluster;
00797 
00798       masterClusterItr.ResetFirst();
00799       
00800       //set the iterator direction to fill the vector of hit planes
00801       if(idir == 1){
00802         clusterItr.ResetFirst();
00803         planeClusterItr.ResetFirst();
00804       }
00805       else if(idir == -1){
00806         clusterItr.ResetLast();
00807         planeClusterItr.ResetLast();
00808       }
00809       
00810       //create a vector of ints to hold the numbers of 
00811       //the planes with acceptable clusters in them
00812       //only holds the numbers for  planes in the current view
00813       vector<Int_t> hitPlanes;
00814       Int_t prevPlane = 0;
00815       while( (cluster = clusterItr()) ){
00816         
00817         MSG("TrackSR", Msg::kDebug) << "cluster on plane " 
00818                                    << cluster->GetPlane() << " time = " 
00819                                    << cluster->GetBegTime() << endl;
00820         if(cluster->GetPlane() != prevPlane){
00821           hitPlanes.push_back(cluster->GetPlane());
00822           prevPlane = cluster->GetPlane();
00823         }
00824         
00825       }//end loop over clusters in the iterator
00826       
00827       FindTimingDirection(clusterItr);
00828       
00829       Int_t iterate = 0;
00830       bool isDone = false;
00831       Int_t clustersLeft = clusterItr.SizeSelect();
00832       Int_t begPlane = 0;
00833       Int_t loopTracks = 1;
00834       
00835       //do the following loop until the track is done or
00836       //there are no more valid clusters left
00837       //     while( !isDone && clustersLeft>0 && loopTracks>0){
00838       while( !isDone && clustersLeft>0){
00839         
00840         //set the iterator direction to fill the vector of hit planes
00841         if(idir == 1){
00842           clusterItr.ResetFirst();
00843           planeClusterItr.ResetFirst();
00844         }
00845         else if(idir == -1){
00846           clusterItr.ResetLast();
00847           planeClusterItr.ResetLast();
00848         }
00849         Track2DSR *track = 0;
00850         Int_t trackCtr = 0;     
00851         //reset the number of clusters left in this view
00852         clustersLeft = 0;
00853         
00854         //identify the seed clusters and put them into trackList
00855         loopTracks = FillTrackList(clusterItr, trackList, seedClusters,
00856                                    houghSlope, houghIntercept, idir, iterate, 
00857                                    begPlane);
00858         
00859         //get an iterator over the tracks
00860         Track2DSRItr trackItr(trackList);
00861         Track2DSRItr trackItr2(trackList);
00862         
00863         //program the selection key functions to get only tracks
00864         //in the current view
00865         Track2DSRKeyFunc *select1KF = trackItr.CreateKeyFunc();
00866         Track2DSRKeyFunc *select2KF = trackItr2.CreateKeyFunc();
00867         
00868         if(viewCtr == 0){
00869           select1KF->SetFun(KeyOnUViewTracks);
00870           select2KF->SetFun(KeyOnUViewTracks);
00871         }
00872         else if(viewCtr == 1){
00873           select1KF->SetFun(KeyOnVViewTracks);
00874           select2KF->SetFun(KeyOnVViewTracks);
00875         }
00876         trackItr.GetSet()->AdoptSelectKeyFunc(select1KF);
00877         trackItr2.GetSet()->AdoptSelectKeyFunc(select2KF);
00878         
00879         select1KF = 0;
00880         select2KF = 0;
00881 
00882 
00883         MSG("TrackSR",Msg::kDebug) << "after FillTrackList, total # tracks = " 
00884                                   << trackList->GetLast()+1 << "\n";
00885 
00886         //loop over the hit planes and tracks made in this iteration
00887         //to add clusters to the tracks
00888         AddClustersToTracks(trackItr, planeClusterItr, idir, iterate, 
00889                             hitPlanes,trk2dmin);
00890 
00891 
00892         MSG("TrackSR",Msg::kDebug) << " after add clusters, total # tracks = " 
00893                                   << trackList->GetLast()+1 << "\n";
00894 
00895         IdentifyBadTracks(trackItr, iterate, trk2dmin);
00896 
00897           trackItr.ResetFirst();
00898           trackCtr = 0;
00899           while( (track = trackItr()) ){
00900 
00901             MSG("TrackSR", Msg::kDebug) << "track " << trackCtr 
00902                                        << " is bad " << (Int_t)track->IsBad() << endl; 
00903             for(Int_t i = 0; i<=track->GetLast(); ++i){
00904               cluster = track->GetCluster(i);
00905               MSG("TrackSR", Msg::kDebug) << "\tTC plane = " << cluster->GetPlane()
00906                                          << " tpos = " << cluster->GetTPos() 
00907                                          << endl;
00908             }
00909             ++trackCtr;
00910           }     
00911         MarkUsedClusters(trackItr, clusterItr, iterate, nmaxtrack, houghSlope, 
00912                          houghIntercept);
00913 
00914  
00915         //subtract the # of tracks from the count of that are formed and then
00916         //removed in this loop
00917         loopTracks -= RemoveTrackSubsets(trackItr, trackItr2, trackList);
00918          
00919         MSG("TrackSR",Msg::kDebug) << "after remove subsets, total # tracks = " 
00920                                   << trackList->GetLast()+1 << "\n";
00921           
00922         //need to count the number of valid clusters left in the clusterItr
00923         clusterItr.ResetFirst();
00924         while( (cluster = clusterItr()) ){
00925           if(cluster->IsValid()){
00926             ++clustersLeft;
00927             MSG("TrackSR", Msg::kDebug) << "valid cluster on plane " 
00928                                        << cluster->GetPlane() << endl;
00929           }
00930         }
00931         
00932         //see if the tracks span the event
00933         // if cosmic mode and found tracks span event stop
00934         if(fParmIsCosmic && sliceVtxPlane>0 
00935            && sliceEndPlane>0 ){
00936           Bool_t isDoneU=false;
00937           Bool_t isDoneV=false;
00938           for(Int_t itrk=0; itrk<=trackList->GetLast(); itrk++){
00939             Track2DSR *track = dynamic_cast<Track2DSR*>(trackList->At(itrk));
00940             
00941             //if the number of clusters in the track is the same as 
00942             //the number of hit planes in the view, stop tracking
00943             if(track->GetLast()+1==sliceNPlane){
00944               if(track->GetPlaneView()==PlaneView::kU) isDoneU = true;
00945               if(track->GetPlaneView()==PlaneView::kV) isDoneV = true;
00946               isDone = isDoneU & isDoneV;
00947               if(isDone)break;
00948             }
00949           }//end loop over tracks
00950         }//end if a cosmic event 
00951 
00952         //another iteration down, increment the counter
00953         ++iterate;
00954       }//end loop to find all the tracks in the view
00955             
00956 
00957      //now fill in the gaps in the tracks for this view
00958       Track2DSRItr trackItr(trackList);
00959       FillInGaps(trackItr, masterClusterItr, idir);
00960       
00961       //clear the iterator select functions for the next view
00962       masterClusterItr.GetSet()->AdoptSelectKeyFunc(0);
00963       planeClusterItr.GetSet()->AdoptSelectKeyFunc(0);
00964       clusterItr.GetSet()->AdoptSelectKeyFunc(0);
00965       
00966       //increment the view cnt to do the next view
00967       ++viewCtr;
00968     }//end loop over each view for 2d tracking
00969     
00970     // merge tracks which are non-overlapping, and which are consistent
00971     MergeTracks(trackList,idir);
00972     
00973     //get iterators over the u and v view tracks
00974     Track2DSRItr uTrackItr(trackList);
00975     Track2DSRItr vTrackItr(trackList);
00976     
00977     MSG("TrackSR", Msg::kDebug) << "number of 2d tracks = " << trackList->GetLast()+1  << endl;
00978     
00979     CandSliceHandle * trackSlice=slice;
00980     
00981     // ********************************************
00982     // for near detector, add spectrometer hits to tracks which end near
00983     // plane 121. 
00984     // *******************************************
00985     if(EnableSpectFind){
00986       if(fDetector==Detector::kNear){
00987         CandSliceHandle *  dupSlice=slice->DupHandle();
00988         SpectrometerTracking(trackList, cth, slice, dupSlice, cx);
00989         // if spectrometer tracking has cloned slice, associate that slice
00990         // with track
00991         if(dupSlice->GetUidInt()!=slice->GetUidInt()){
00992           trackSlice = dynamic_cast<CandSliceHandle *>
00993             (fCSlcLHnew->FindDaughter(dupSlice));
00994         }
00995         delete dupSlice;
00996       }  
00997     }
00998     //make the 3d tracks out of the 2d tracks
00999     TObjArray *candTracks = new TObjArray(1,0);
01000  
01001     //program the selection key functions
01002     Track2DSRKeyFunc *selectUKF = uTrackItr.CreateKeyFunc();
01003     Track2DSRKeyFunc *selectVKF = vTrackItr.CreateKeyFunc();
01004 
01005     
01006     selectUKF->SetFun(KeyOnUViewTracks);
01007     selectVKF->SetFun(KeyOnVViewTracks);
01008     
01009     uTrackItr.GetSet()->AdoptSelectKeyFunc(selectUKF);
01010     vTrackItr.GetSet()->AdoptSelectKeyFunc(selectVKF);
01011     
01012     selectUKF = 0;
01013     selectVKF = 0;
01014   
01015     FormCandTrackSR(uTrackItr, vTrackItr, candTracks, ch, trackSlice, ah, cxx, idir);
01016     
01017     //get iterators over the cand tracks
01018     CandTrackSRHandleItr candTracks1(candTracks);
01019     CandTrackSRHandleItr candTracks2(candTracks);
01020     
01021     //get rid of possible twin tracks
01022     DeleteTwinTracks(candTracks1, candTracks2, candTracks,
01023                      maxUTrack, maxVTrack,nmax3dtrk, ch);
01024     
01025 
01026     //free up the memory used by trackList and trackClusterList
01027     CleanUp(trackList, trackClusterList);    
01028     delete trackClusterList;
01029     delete trackList;
01030     delete seedClusters;
01031     delete candTracks;
01032     
01033   }//end loop over slices passed in
01034 
01035   // remove strips from striplist in spectrometer
01036   // if not on any track
01037 
01038   if(EnableSpectFind){
01039     if(fDetector==Detector::kNear){
01040       RemoveUnusedSpectStrips(cth,fCSLHnew);
01041       RemoveStripsInSlice(cth,fCSlcLHnew);
01042     }
01043   }
01044 
01045   if(fCDLHnew){
01046     CandRecord *crec = cx.GetCandRecord();
01047     assert(crec);
01048     fCDLHnew->SetName("canddigitlist");
01049     crec->SecureCandHandle(*fCDLHnew);
01050     delete fCDLHnew;
01051   }
01052   if(fCSlcLHnew){
01053     CandRecord *crec = cx.GetCandRecord();
01054     assert(crec);
01055     fCSlcLHnew->SetName("CandSliceList");
01056     crec->SecureCandHandle(*fCSlcLHnew);
01057     delete fCSlcLHnew;
01058   }
01059   if(fCSLHnew){ 
01060     CandRecord *crec = cx.GetCandRecord();
01061     assert(crec);
01062     fCSLHnew->SetName("CandStripList");
01063     crec->SecureCandHandle(*fCSLHnew);
01064     delete fCSLHnew;
01065   }
01066 
01067   MSG("TrackSR", Msg::kDebug) << "done with AlgTrackSRList" << endl;
01068 
01069   return;
01070 }//end Run Alg

void AlgTrackSRList::SpectrometerTracking ( TObjArray *  trackList,
CandTrackSRListHandle cth,
CandSliceHandle slice,
CandSliceHandle dupSlice,
CandContext  cx 
) [private]

Definition at line 1178 of file AlgTrackSRList.cxx.

References Track2DSR::Add(), CandHandle::AddDaughterLink(), TrackClusterSR::AddStrip(), CandTrackSRListHandle::AddTrackCluster(), Track2DSR::Compress(), CandDigitListHandle::DupHandle(), CandStripListHandle::DupHandle(), CandSliceListHandle::DupHandle(), CandDigitHandle::DupHandle(), fCDLHnew, fCSlcLHnew, fCSLHnew, fHoughUIntercept, fHoughUSlope, fHoughVIntercept, fHoughVSlope, CandRecord::FindCandHandle(), CandHandle::FindDaughter(), FindNumSkippedPlanes(), fParmIsCosmic, fParmMisalignmentError, fParmTrk2DNContiguous, fParmTrk2DNSkip, fParmTrk2DNSkipHit, fParmTrk2DNSkipRemove, fParmTrkClsNSkip, fSliceEndPlane, fSliceEndUPlane, fSliceEndVPlane, Track2DSR::GetBegPlane(), PlexSEIdAltL::GetBestItem(), PlexSEIdAltL::GetBestWeight(), CandHandle::GetCandRecord(), CandContext::GetCandRecord(), Track2DSR::GetCluster(), CandHandle::GetDaughterIterator(), Track2DSR::GetEndPlane(), Track2DSR::GetForwardSlope(), MomNavigator::GetFragment(), AlgFactory::GetInstance(), Track2DSR::GetLast(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), CandContext::GetMom(), TrackClusterSR::GetNStrip(), TrackClusterSR::GetPlane(), VHS::GetPlane(), CandStripHandle::GetPlane(), TrackClusterSR::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), CandDigitHandle::GetPlexSEIdAltLWritable(), PlexSEIdAltLItem::GetSEId(), VHS::GetStrip(), CandStripHandle::GetStrip(), PlexStripEndId::GetStrip(), TrackClusterSR::GetStripList(), TrackClusterSR::GetTPos(), CandHandle::GetUidInt(), RecMinos::GetVldContext(), TrackClusterSR::GetZPos(), IsBestClusterInPlane(), TrackClusterSR::IsDoubleEnded(), Msg::kDebug, TrackClusterSR::KeyFromPlane(), PlaneView::kU, PlaneView::kUnknown, PlaneView::kV, CandStrip::MakeCandidate(), MSG, Track2DSR::RemoveAt(), CandHandle::RemoveDaughter(), CandHandle::SetCandRecord(), CandHandle::SetName(), and tc.

Referenced by RunAlg().

01183 {
01184   MSG("SpectTrackSR",Msg::kDebug) << " In Spectrometer Tracking " << endl;
01185   // Create Candcontext
01186   CandContext cxx(this,cx.GetMom());
01187   
01188   // Get singleton instance of AlgFactory
01189   AlgFactory &af = AlgFactory::GetInstance();
01190   
01191   CandRecord *candrec = cx.GetCandRecord();
01192   assert(candrec);
01193   const VldContext *vldcptr = candrec->GetVldContext();
01194   assert(vldcptr);
01195   VldContext vldc = *vldcptr;
01196   UgliGeomHandle ugh(vldc);
01197 
01198   // generate clones of the digitlist and striplist
01199   const MomNavigator * mom = cx.GetMom();
01200   CandRecord *crec = dynamic_cast<CandRecord *>
01201     (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
01202 
01203   CandStripListHandle *cslh = dynamic_cast<CandStripListHandle *>
01204     (crec->FindCandHandle("CandStripListHandle"));
01205   assert(cslh);
01206   CandSliceListHandle *cslclh = dynamic_cast<CandSliceListHandle *>
01207     (crec->FindCandHandle("CandSliceListHandle"));
01208   assert(cslclh);
01209   CandDigitListHandle *cdlh = dynamic_cast<CandDigitListHandle *>
01210     (crec->FindCandHandle("CandDigitListHandle","canddigitlist"));
01211   assert(cdlh);
01212   
01213   if(!fCDLHnew){
01214     fCDLHnew = cdlh->DupHandle();
01215     fCDLHnew->SetCandRecord(cdlh->GetCandRecord());
01216     fCDLHnew->SetName("candspecdigitlist");
01217   }
01218   if(!fCSlcLHnew){
01219     fCSlcLHnew = cslclh->DupHandle();
01220     fCSlcLHnew->SetCandRecord(cslclh->GetCandRecord());
01221     fCSlcLHnew->SetName("candspecslicelist");
01222   }
01223   if(!fCSLHnew){
01224     fCSLHnew = cslh->DupHandle();     
01225     fCSLHnew->SetCandRecord(cslh->GetCandRecord());
01226     fCSLHnew->SetName("candspecstriplist");
01227   }
01228 
01229   TIter newdigitItr(fCDLHnew->GetDaughterIterator());
01230   CandStripHandleItr oldstripItr(cslh->GetDaughterIterator());
01231   
01232   Int_t specidir=1;
01233   //  Int_t minview = min(PlaneView::kU,PlaneView::kV);
01234   
01235   oldstripItr.Reset();
01236   oldstripItr.SetDirection(specidir);
01237   CandStripHandleItr  newstripItr(fCSLHnew->GetDaughterIterator());     
01238 
01239 // generate cluster list for spectrometer hits
01240 // clusters are made for all possible strip end alternatives
01241 
01242 //this map stores the correspondence between the new strips we will
01243 // create here and the originals.
01244 
01245   map<CandStripHandle*,CandStripHandle*>newstripmap;
01246   map<CandStripHandle*,CandStripHandle*>slicemap;
01247 
01248   Int_t specbegplane=1000; 
01249   TObjArray spectrackclusterlist;
01250 
01251   AlgHandle stripah = af.GetAlgHandle("AlgStripSR","default");
01252        
01253   CandStripHandleItr slicestripItr(slice->GetDaughterIterator());
01254   for (CandStripHandle *strip = slicestripItr(); strip ; strip = slicestripItr()) {
01255     if (strip->GetPlane()>120 && (strip->GetPlaneView()==PlaneView::kU || strip->GetPlaneView()==PlaneView::kV) ) {
01256       MSG("SpectTrackSR",Msg::kDebug) << " spect strip " << strip->GetPlane() << " " << strip->GetStrip() << " " << strip->GetCharge() << endl;
01257 
01258       if(strip->GetPlaneView()==PlaneView::kU && strip->GetPlane()>fSliceEndUPlane)
01259         fSliceEndUPlane=strip->GetPlane();
01260       if(strip->GetPlaneView()==PlaneView::kV && strip->GetPlane()>fSliceEndVPlane)
01261         fSliceEndVPlane=strip->GetPlane();
01262 
01263       TIter digitItr(strip->GetDaughterIterator());
01264       CandDigitHandle *dig=
01265         dynamic_cast<CandDigitHandle*>(digitItr());           
01266       const PlexSEIdAltL& altl = dig->GetPlexSEIdAltL();
01267       int siz = altl.size();
01268       for (int ind = 0; ind < siz; ++ind) {
01269         TObjArray diglist;  
01270         TIter newstripdauItr(strip->GetDaughterIterator());
01271         while( CandDigitHandle * olddig = dynamic_cast<CandDigitHandle*>(newstripdauItr())){
01272           CandDigitHandle * newdig = olddig->DupHandle();
01273           PlexSEIdAltL& newaltl = newdig->GetPlexSEIdAltLWritable(); 
01274           for (int newind = 0; newind < siz; ++newind) {
01275             if(newind==ind){
01276               newaltl[newind].SetWeight((float)1.);
01277             }
01278             else{
01279               newaltl[newind].SetWeight((float)0.);
01280             }
01281           }
01282           newdig->SetCandRecord(olddig->GetCandRecord());
01283           diglist.Add(newdig);   // diglist does not own newdig.  These must be explicitely deleted 
01284         }
01285         cxx.SetDataIn(&diglist);
01286         CandStripHandle newstrip =
01287           CandStrip::MakeCandidate(stripah,cxx);
01288         newstrip.SetCandRecord(cslh->GetCandRecord());
01289         CandStripHandle *daustrip = new CandStripHandle(newstrip);
01290         CandStripHandle * nwstrip = dynamic_cast<CandStripHandle *>
01291           (fCSLHnew->FindDaughter(strip));
01292 
01293         for (Int_t i=0; i<=diglist.GetLast(); i++) {
01294           CandDigitHandle * deldig =  dynamic_cast<CandDigitHandle*>(diglist.At(i));
01295           delete deldig;
01296         }
01297         
01298         Bool_t found=false;
01299         Int_t listlen = spectrackclusterlist.GetLast();
01300         for (Int_t i=0; i<=listlen; i++) {
01301           TrackClusterSR* oldcluster = dynamic_cast<TrackClusterSR*>(spectrackclusterlist.At(i));
01302           assert(oldcluster);
01303           Int_t minstrip = oldcluster->GetMinStrip();
01304           Int_t maxstrip = oldcluster->GetMaxStrip();
01305           if (newstrip.GetPlane()==oldcluster->GetPlane() &&
01306               (TMath::Abs(minstrip-newstrip.GetStrip())<=fParmTrkClsNSkip+1 || 
01307                TMath::Abs(maxstrip-newstrip.GetStrip())<=fParmTrkClsNSkip+1)){
01308             oldcluster->AddStrip(daustrip);
01309             Int_t istrip=oldcluster->GetNStrip()-1; 
01310             CandStripHandle * clusterstrip = dynamic_cast<CandStripHandle *>(oldcluster->GetStripList()->At(istrip));
01311             slicemap[clusterstrip]=strip;
01312             newstripmap[clusterstrip]=nwstrip;
01313             found=true;
01314             break;
01315           }
01316         }
01317         if(!found) {
01318           MSG("SpectTrackSR",Msg::kDebug) << " spect cluster " << newstrip.GetPlane() << endl;    
01319           TrackClusterSR *trackcluster = new TrackClusterSR(daustrip, 
01320                                                             fParmMisalignmentError);
01321           Int_t istrip=trackcluster->GetNStrip()-1; 
01322           CandStripHandle * clusterstrip = dynamic_cast<CandStripHandle *>(trackcluster->GetStripList()->At(istrip));
01323           slicemap[clusterstrip]=strip;
01324           newstripmap[clusterstrip]=nwstrip;
01325           spectrackclusterlist.Add(trackcluster);
01326         
01327           if (newstrip.GetPlane()<specbegplane)specbegplane = 
01328                                                        newstrip.GetPlane();
01329         }
01330         delete daustrip;
01331       }
01332     }
01333   }   
01334   fSliceEndPlane = TMath::Max(fSliceEndUPlane,fSliceEndVPlane);
01335 
01336   // now remove all CandStrips in the spectrometer in the slice list
01337   // and the new strip list.  
01338        
01339   slicestripItr.Reset();
01340   
01341   MSG("SpectTrackSR",Msg::kDebug) << "# spectrometer track clusters: "
01342                              << spectrackclusterlist.GetLast() << endl;
01343   MSG("SpectTrackSR",Msg::kDebug) << "eliminating low pulse height spectrometer clusters" 
01344                              << endl;
01345   
01346  
01347   TrackClusterSRItr spectclusterItr(&spectrackclusterlist);
01348   TrackClusterSRKeyFunc *spectclusterKf = spectclusterItr.CreateKeyFunc();
01349   spectclusterKf->SetFun(TrackClusterSR::KeyFromPlane);
01350   spectclusterItr.GetSet()->AdoptSortKeyFunc(spectclusterKf);
01351   spectclusterKf = 0;
01352   spectclusterItr.SetDirection(specidir);
01353   spectclusterItr.Reset();
01354   
01355   // create vector of hit planes in the spectrometer
01356   
01357   vector<Int_t> specplaneindex;
01358   vector<PlaneView::PlaneView_t> specplaneviewindex;
01359   Int_t ncount=0;
01360   Int_t oldipln=0;
01361   while (TrackClusterSR *tc = spectclusterItr()) {
01362     Int_t ipln = tc->GetPlane();
01363     if (!ncount || ipln!=oldipln) {
01364       specplaneindex.push_back(ipln);
01365       specplaneviewindex.push_back(tc->GetPlaneView());
01366     }
01367     ncount++;
01368     oldipln = ipln;
01369   }
01370   
01371   // now  add spectrometer track hits to existing tracks
01372   // no new tracks are formed at this point
01373   
01374   Bool_t trackcoverageu(1);
01375   Bool_t trackcoveragev(1);
01376   map<Track2DSR*,Int_t> nhitplaneskip;
01377   
01378   // loop over hit planes in the spectrometer, from front to back
01379   
01380   for (Int_t iplnindx=0; iplnindx<static_cast<Int_t>(specplaneindex.size()) && (trackcoverageu || trackcoveragev);iplnindx++) {
01381     Int_t ipln = specplaneindex[iplnindx];
01382     MSG("SpectTrackSR",Msg::kDebug) << "PLANE " << ipln << "/" << specbegplane << "\n";
01383     if (ipln>=specbegplane) {
01384       
01385       Bool_t hasdoubleended(0);
01386       
01387       // loop over spectrometer clusters, finding those in plane ipln
01388       // these clusters are held in tclusterpln
01389       
01390       TObjArray tclusterpln;
01391       for (Int_t i=0; i<=spectrackclusterlist.GetLast(); i++) {
01392         TrackClusterSR* tc = dynamic_cast<TrackClusterSR*>(spectrackclusterlist.At(i));
01393         if (!hasdoubleended && tc->IsDoubleEnded()) {
01394           hasdoubleended = 1;
01395         }
01396         
01397         if (tc->GetPlane()==ipln) {
01398           tclusterpln.Add(tc);
01399           if (tc->GetPlaneView()==PlaneView::kU) {
01400             trackcoverageu = 0;
01401           }
01402           if (tc->GetPlaneView()==PlaneView::kV) {
01403             trackcoveragev = 0;
01404           }
01405         }
01406       }
01407       Int_t useph = 0;
01408       if (fParmIsCosmic==1 && hasdoubleended) {
01409         useph = 1;
01410       }
01411       // loop over tracks. Each track will be extrapolated to and compared
01412       // with the position of each cluster on plane ipln
01413       
01414       for (Int_t itrk=0; itrk<=track2dlist->GetLast(); itrk++) {
01415         Track2DSR *thistrack = dynamic_cast<Track2DSR*>(track2dlist->At(itrk));
01416         MSG("SpectTrackSR",Msg::kDebug) << "track " << itrk << "/" << track2dlist->GetLast() << "\n";
01417         Double_t residParams[] = {1.e6,1.e6};
01418         TrackClusterSR *besttc=0;
01419         PlaneView::PlaneView_t tcplaneview = PlaneView::kUnknown;
01420         Bool_t shouldhit(kFALSE); // true is same view as track
01421         
01422         // loop over track end interval to extrapolate from
01423         //  Int_t ibefpln=0;
01424         for (Int_t ibefpln=0; ibefpln<=fParmTrk2DNSkipRemove && ibefpln<thistrack->GetLast() && !besttc; ibefpln++) {
01425           Int_t numSkippedHit = 0;
01426           TrackClusterSR *tct = 
01427             (thistrack->GetCluster(ibefpln));
01428           
01429           Int_t slopelookup=ibefpln;
01430           if(slopelookup<thistrack->GetLast())slopelookup++;
01431           Double_t trackSlope = thistrack->GetForwardSlope(ibefpln+1);
01432           Int_t deltaPlane=0;
01433           Int_t dview=0;
01434           Int_t numSkippedActive=0; // number of active planes of same view
01435           
01436           // loop over all clusters in plane ipln 
01437           
01438           for (Int_t itc=0; itc<=tclusterpln.GetLast(); itc++) {
01439             TrackClusterSR *tc = dynamic_cast<TrackClusterSR*>(tclusterpln.At(itc));
01440             
01441             if(!itc){
01442               deltaPlane = tc->GetPlane()-tct->GetPlane();
01443               tcplaneview = tc->GetPlaneView();
01444               dview = tcplaneview-tct->GetPlaneView();
01445               
01446               if (dview==0) shouldhit = kTRUE;
01447               if (shouldhit && deltaPlane>0) {
01448                 Float_t upos=0;
01449                 Float_t vpos=0;
01450                 Float_t uslope=0;
01451                 Float_t vslope=0;
01452                 if(tcplaneview==PlaneView::kU){
01453                   upos = tct->GetTPos();
01454                   vpos = fHoughVIntercept + tct->GetZPos()*fHoughVSlope;
01455                   uslope=trackSlope;
01456                   vslope=fHoughVSlope;
01457                 }
01458                 else{
01459                   vpos = tct->GetTPos();
01460                   upos = fHoughUIntercept + tct->GetZPos()*fHoughUSlope;
01461                   vslope=trackSlope;
01462                   uslope=fHoughUSlope;
01463                 }
01464                 numSkippedActive=FindNumSkippedPlanes(tc->GetPlane(),
01465                                                       tct->GetPlane(),
01466                                                       upos,vpos,
01467                                                       tct->GetZPos(),
01468                                                       uslope,
01469                                                       vslope,
01470                                                       specidir,
01471                                                       tct->GetPlaneView());
01472                 
01473               }
01474             }
01475             Int_t nhit = thistrack->GetLast()-ibefpln+1;
01476             Float_t aveSpacing=2;
01477             if(thistrack->GetLast()>0){
01478               aveSpacing = TMath::Abs(thistrack->GetEndPlane()-thistrack->GetBegPlane())/ (thistrack->GetLast());
01479             }         
01480             MSG("SpectTrackSR",Msg::kDebug) << "nhit, nplaneskip = " << nhit << "/" << fParmTrk2DNContiguous  
01481                                             << "  " << numSkippedActive << "/" << fParmTrk2DNSkip 
01482                                             << " dview/dplane " << dview << "/" << deltaPlane 
01483                                             << " track plane/min track plane " << tct->GetPlane() 
01484                                             << "/" << 120 - aveSpacing*fParmTrk2DNSkip*2 << endl;
01485             
01486             if (dview==0 && 
01487                 deltaPlane>=0 &&  
01488                 numSkippedActive<=fParmTrk2DNSkip && 
01489                 //  tct->GetPlane()>=(120-aveSpacing*fParmTrk2DNSkip*2) && 
01490                 (nhit>=fParmTrk2DNContiguous || numSkippedActive==0)) { 
01491               
01492               MSG("SpectTrackSR",Msg::kDebug) << " checking SE alternative: Tpos: " << tc->GetTPos() <<  endl;
01493               if (tcplaneview==PlaneView::kU) {
01494                 trackcoverageu = 1;
01495                 MSG("SpectTrackSR",Msg::kDebug) << "clearing U trackcoverage\n";
01496               }
01497               if (tcplaneview==PlaneView::kV) {
01498                 trackcoveragev = 1;
01499                 MSG("SpectTrackSR",Msg::kDebug) << "clearing V trackcoverage\n";
01500               }
01501               
01502               //reset the residParams because you are on a new plane
01503               Bool_t isSpectrometerPlane=true;
01504               
01505               if(IsBestClusterInPlane(tc,tct, numSkippedHit,
01506                                       nhit, ibefpln, trackSlope,thistrack, 
01507                                       
01508                                       residParams, 
01509                                       deltaPlane,isSpectrometerPlane))  besttc = tc;
01510             }
01511           }
01512                 } // end ibefpln
01513         
01514         if (besttc) { 
01515           TrackClusterSR * newtrackcluster = 0;
01516           
01517           // loop over strips in this cluster, and for each....
01518           //    loop over digit daughters, and find on new digit list
01519           //     3 possibilities exist...
01520           //     new digit has zero weights - clone with correct weights
01521           //         and delete zero weight version.
01522           //     new digit has weights already set (correctly) do nothing
01523           //     new digit has weight set, but pointing to a different strip
01524           //         in this case, make new CandDigit.  
01525           
01526           for (Int_t istrip=0; istrip<besttc->GetNStrip(); istrip++) {
01527             CandStripHandle * strip = dynamic_cast<CandStripHandle *>(besttc->GetStripList()->At(istrip));
01528             assert(strip);
01529             
01530             // now build a new strip to be added to strip list
01531             
01532             TObjArray cdhAr;
01533             cdhAr.Clear();
01534             Bool_t deleteold=true;
01535             // loop over digit daughters of strip
01536             
01537             TIter digitItr(strip->GetDaughterIterator());
01538             while (CandDigitHandle *dig=dynamic_cast<CandDigitHandle*>(digitItr())) {    
01539               newdigitItr.Reset();
01540               Bool_t moddig=false;
01541               Bool_t prevadded=false;
01542               CandDigitHandle * founddig=0;
01543               
01544               // find this digit in the new digit list
01545               
01546               while (CandDigitHandle *newdig=dynamic_cast<CandDigitHandle*>(newdigitItr())){
01547                 if(dig->GetChannelId()==newdig->GetChannelId() && dig->GetRawDigitIndex()==newdig->GetRawDigitIndex()){ 
01548                   founddig=newdig;                     
01549                   const PlexSEIdAltL& oldaltl = newdig->GetPlexSEIdAltL();
01550                   
01551                   // digit is found - check weights
01552                   
01553                   if(oldaltl.GetBestWeight()==0) {
01554                     // if best weight is zero, we make a version of this
01555                     // digit with correct weights first duplicating the handle.
01556                     CandDigitHandle * dupdig = founddig->DupHandle(); 
01557                     PlexSEIdAltL& dupaltl = dupdig->GetPlexSEIdAltLWritable();
01558                     moddig=true;
01559                     int siz = dupaltl.size();
01560                     for (int ind = 0; ind < siz; ++ind) {
01561                       if(dupaltl[ind].GetSEId().GetStrip()==
01562                          strip->GetStrip()){
01563                         dupaltl[ind].SetWeight((float)1.);
01564                       }
01565                       else{
01566                         dupaltl[ind].SetWeight((float)0.);
01567                       }
01568                     }
01569                     
01570                     dupdig->SetCandRecord(founddig->GetCandRecord());
01571                     fCDLHnew->AddDaughterLink(*dupdig);
01572                     CandDigitHandle * daudig = dynamic_cast<CandDigitHandle *>(fCDLHnew->FindDaughter(dupdig));
01573                     assert(daudig);
01574                     fCDLHnew->RemoveDaughter(newdig);
01575                     cdhAr.Add(daudig);
01576                     delete dupdig;
01577                   }
01578                   else {
01579                     
01580                     // this digit has already been set to the correct SEId, 
01581                     // so just leave it alone.
01582                     
01583                     if(oldaltl.GetBestWeight()==1.0 && oldaltl.GetBestItem().GetSEId().GetStrip()==strip->GetStrip()){
01584                       prevadded=true;
01585                       cdhAr.Add(newdig);
01586                     }
01587                   }
01588                   break;
01589                 }
01590               }
01591               
01592               if(founddig && !moddig && !prevadded){ 
01593                 // need to construct new CandDigit, as there is
01594                 // no existing digit which has proper altl weights.  Since
01595                 // the weights have apparently already been set in the
01596                 // found digit (probably associated 
01597                 // with an earlier track),
01598                 // leave it in place.
01599                 deleteold=false;
01600                 CandDigitHandle * dupdig = founddig->DupHandle(); 
01601                 PlexSEIdAltL& dupaltl = dupdig->GetPlexSEIdAltLWritable();
01602                 int siz = dupaltl.size();
01603                 for (int ind = 0; ind < siz; ++ind) {
01604                   if(dupaltl[ind].GetSEId().GetStrip()==
01605                      strip->GetStrip()){
01606                     dupaltl[ind].SetWeight((float)1.);
01607                   }
01608                   else{
01609                     dupaltl[ind].SetWeight((float)0.);
01610                   }
01611                 }
01612                 dupdig->SetCandRecord(founddig->GetCandRecord());
01613                 fCDLHnew->AddDaughterLink(*dupdig);
01614                 CandDigitHandle * daudig = dynamic_cast<CandDigitHandle *>(fCDLHnew->FindDaughter(dupdig));
01615                 if(daudig)cdhAr.Add(daudig);
01616                 delete dupdig;
01617               }       
01618             }
01619             
01620             // we now construct the strip, and add to striplist and slice
01621             // daughter list
01622             
01623             AlgHandle stripah = af.GetAlgHandle("AlgStripSR","default");
01624             if(cdhAr.GetLast()>=0){
01625               cxx.SetDataIn(&cdhAr);
01626               CandStripHandle striphandle =
01627                 CandStrip::MakeCandidate(stripah,cxx);
01628               striphandle.SetCandRecord(cslh->GetCandRecord());
01629 
01630               fCSLHnew->AddDaughterLink(striphandle);           
01631               dupSlice->AddDaughterLink(striphandle);
01632               CandStripHandle * daustrip = dynamic_cast<CandStripHandle *>(fCSLHnew->FindDaughter(&striphandle));
01633 
01634               assert (daustrip);
01635               // add to track cluster list
01636               if(istrip==0 || !newtrackcluster){
01637                 newtrackcluster = new TrackClusterSR(daustrip, 
01638                                                      fParmMisalignmentError);
01639               }
01640               else{
01641                 newtrackcluster->AddStrip(daustrip);
01642               }
01643 
01644             }
01645           }
01646           if(newtrackcluster){
01647             // add this cluster to track
01648             TrackClusterSR *tc0 = thistrack->GetCluster(0);
01649             Double_t slope = (newtrackcluster->GetTPos()-tc0->GetTPos())/
01650               (newtrackcluster->GetZPos()-tc0->GetZPos());
01651             // add the full spectrometer cluster to the track, and to the
01652             // track cluster list
01653             thistrack->Add(newtrackcluster,slope);
01654             cth.AddTrackCluster(newtrackcluster);
01655             delete newtrackcluster;
01656           }
01657           nhitplaneskip[thistrack]=0;
01658         }
01659         else if (shouldhit && ipln>thistrack->GetCluster(thistrack->GetLast())->GetPlane()) {
01660           nhitplaneskip[thistrack]++;
01661         }
01662       }
01663     }
01664   }
01665   spectrackclusterlist.Delete();
01666 
01667   if(dupSlice->GetUidInt()!=slice->GetUidInt()){
01668     fCSlcLHnew->AddDaughterLink(*dupSlice);
01669     fCSlcLHnew->RemoveDaughter(slice);
01670   }
01671   // check each track for outlier hit at two track ends, and remove if found
01672   for (Int_t itrk=0; itrk<=track2dlist->GetLast(); itrk++) {
01673     Track2DSR *track = dynamic_cast<Track2DSR*>(track2dlist->At(itrk)); 
01674     if(track->GetLast()>2){
01675       Float_t gap1=TMath::Abs(track->GetCluster(0)->GetPlane()-track->GetCluster(1)->GetPlane());
01676       Float_t prevgap1=TMath::Abs(track->GetCluster(2)->GetPlane()-track->GetCluster(1)->GetPlane());
01677       gap1=gap1/prevgap1;
01678       Float_t gap2=TMath::Abs(track->GetCluster(track->GetLast())->GetPlane()-track->GetCluster(track->GetLast()-1)->GetPlane());
01679       Float_t prevgap2=TMath::Abs(track->GetCluster(track->GetLast()-1)->GetPlane()-track->GetCluster(track->GetLast()-2)->GetPlane());
01680       gap2=gap2/prevgap2;
01681       Int_t gapfactor=2;
01682       if(gap1>fParmTrk2DNSkipHit*gapfactor)track->RemoveAt(0);
01683       if(gap2>fParmTrk2DNSkipHit*gapfactor)track->RemoveAt(track->GetLast());
01684       track->Compress();
01685     }   
01686   }
01687 }

void AlgTrackSRList::Trace ( const char *  c  )  const [virtual]

Reimplemented from AlgBase.

Definition at line 4449 of file AlgTrackSRList.cxx.

04450 {
04451 }


Member Data Documentation

Definition at line 192 of file AlgTrackSRList.h.

Referenced by RunAlg(), and SpectrometerTracking().

Definition at line 194 of file AlgTrackSRList.h.

Referenced by RunAlg(), and SpectrometerTracking().

Definition at line 193 of file AlgTrackSRList.h.

Referenced by RunAlg(), and SpectrometerTracking().

Definition at line 198 of file AlgTrackSRList.h.

Referenced by AddClustersToTracks(), MergeTracks(), RunAlg(), and SpectrometerTracking().

Float_t AlgTrackSRList::fHoughUSlope [private]

Definition at line 196 of file AlgTrackSRList.h.

Referenced by AddClustersToTracks(), MergeTracks(), RunAlg(), and SpectrometerTracking().

Definition at line 199 of file AlgTrackSRList.h.

Referenced by AddClustersToTracks(), MergeTracks(), RunAlg(), and SpectrometerTracking().

Float_t AlgTrackSRList::fHoughVSlope [private]

Definition at line 197 of file AlgTrackSRList.h.

Referenced by AddClustersToTracks(), MergeTracks(), RunAlg(), and SpectrometerTracking().

Int_t AlgTrackSRList::fMapIsWide[1000] [private]

Definition at line 119 of file AlgTrackSRList.h.

Referenced by CleanUp(), FindNumSkippedPlanes(), and RunAlg().

Definition at line 169 of file AlgTrackSRList.h.

Referenced by FormCandTrackSR(), and RunAlg().

Definition at line 170 of file AlgTrackSRList.h.

Referenced by FormCandTrackSR(), and RunAlg().

Definition at line 171 of file AlgTrackSRList.h.

Referenced by FormCandTrackSR(), and RunAlg().

Definition at line 172 of file AlgTrackSRList.h.

Referenced by FormCandTrackSR(), and RunAlg().

Definition at line 191 of file AlgTrackSRList.h.

Referenced by FindNumSkippedPlanes(), and RunAlg().

Double_t AlgTrackSRList::fParmHitNTime [private]

Definition at line 137 of file AlgTrackSRList.h.

Referenced by FillTrackList(), IsBestClusterInPlane(), and RunAlg().

Double_t AlgTrackSRList::fParmHitTime [private]

Definition at line 138 of file AlgTrackSRList.h.

Referenced by FillTrackList(), IsBestClusterInPlane(), and RunAlg().

Definition at line 187 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 188 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 189 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 190 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 178 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), RunAlg(), and SpectrometerTracking().

Definition at line 134 of file AlgTrackSRList.h.

Referenced by MakeTrackClusters(), and RunAlg().

Definition at line 133 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 183 of file AlgTrackSRList.h.

Referenced by FindTimingDirection(), and RunAlg().

Definition at line 180 of file AlgTrackSRList.h.

Referenced by MakeTrackClusters(), and RunAlg().

Definition at line 150 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 179 of file AlgTrackSRList.h.

Referenced by Find2DTrackEndPoints(), and RunAlg().

Definition at line 184 of file AlgTrackSRList.h.

Referenced by MakeTrackClusters(), RunAlg(), and SpectrometerTracking().

Definition at line 151 of file AlgTrackSRList.h.

Referenced by FillTrackList(), IsBestClusterInPlane(), and RunAlg().

Double_t AlgTrackSRList::fParmTrk2DAlpha [private]

Definition at line 141 of file AlgTrackSRList.h.

Referenced by AddBestPlaneClusterToTrack(), and RunAlg().

Definition at line 146 of file AlgTrackSRList.h.

Referenced by MergeTracks(), and RunAlg().

Definition at line 145 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 147 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 163 of file AlgTrackSRList.h.

Referenced by CheckForBadClusters(), IdentifyBadTracks(), and RunAlg().

Double_t AlgTrackSRList::fParmTrk2DLinA [private]

Definition at line 158 of file AlgTrackSRList.h.

Referenced by MergeTracks(), and RunAlg().

Double_t AlgTrackSRList::fParmTrk2DLinA0 [private]

Definition at line 154 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Double_t AlgTrackSRList::fParmTrk2DLinB [private]

Definition at line 159 of file AlgTrackSRList.h.

Referenced by MergeTracks(), and RunAlg().

Double_t AlgTrackSRList::fParmTrk2DLinB0 [private]

Definition at line 155 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 142 of file AlgTrackSRList.h.

Referenced by CheckForBadClusters(), FillInGaps(), IdentifyBadTracks(), and RunAlg().

Definition at line 162 of file AlgTrackSRList.h.

Referenced by AddClustersToTracks(), RunAlg(), and SpectrometerTracking().

Definition at line 157 of file AlgTrackSRList.h.

Referenced by CheckForBadClusters(), IsBestClusterInPlane(), and RunAlg().

Definition at line 153 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 161 of file AlgTrackSRList.h.

Referenced by MergeTracks(), and RunAlg().

Definition at line 160 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 143 of file AlgTrackSRList.h.

Referenced by AddClustersToTracks(), RunAlg(), and SpectrometerTracking().

Definition at line 135 of file AlgTrackSRList.h.

Referenced by FillTrackList(), and RunAlg().

Definition at line 167 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 166 of file AlgTrackSRList.h.

Referenced by RemoveTrackSubsets(), and RunAlg().

Definition at line 164 of file AlgTrackSRList.h.

Referenced by RemoveTrackSubsets(), and RunAlg().

Double_t AlgTrackSRList::fParmTrk2DWin0 [private]

Definition at line 136 of file AlgTrackSRList.h.

Referenced by FillTrackList(), and RunAlg().

Definition at line 177 of file AlgTrackSRList.h.

Referenced by DeleteTwinTracks(), and RunAlg().

Definition at line 174 of file AlgTrackSRList.h.

Referenced by DeleteTwinTracks(), and RunAlg().

Definition at line 176 of file AlgTrackSRList.h.

Referenced by MakeTrackClusters(), RunAlg(), and SpectrometerTracking().

Definition at line 144 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 182 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 181 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 131 of file AlgTrackSRList.h.

Referenced by PlaneIsActive().

Definition at line 118 of file AlgTrackSRList.h.

Referenced by FindTimingDirection(), and RunAlg().

Double_t AlgTrackSRList::fSliceDzDs [private]

Definition at line 129 of file AlgTrackSRList.h.

Referenced by IdentifyBadTracks(), and RunAlg().

Definition at line 125 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), RunAlg(), and SpectrometerTracking().

Definition at line 122 of file AlgTrackSRList.h.

Referenced by RunAlg(), and SpectrometerTracking().

Definition at line 123 of file AlgTrackSRList.h.

Referenced by RunAlg(), and SpectrometerTracking().

Definition at line 128 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 126 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 127 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 124 of file AlgTrackSRList.h.

Referenced by IsBestClusterInPlane(), and RunAlg().

Definition at line 120 of file AlgTrackSRList.h.

Referenced by RunAlg().

Definition at line 121 of file AlgTrackSRList.h.

Referenced by RunAlg().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1