AlgTrackSRList.cxx

Go to the documentation of this file.

// $Id: AlgTrackSRList.cxx,v 1.109 2007/02/04 06:03:24 rhatcher Exp $
//
// AlgTrackSRList.cxx
//
// Begin_Html<img src="../../pedestrians.gif" align=center>
// <a href="../source_warning.html">Warning for beginners</a>.<br> 
//
// AlgTrackSRList is a concrete TrackSRList Algorithm class.
//
// Author:  R. Lee 2001.02.26
// Revised: B. Rebel 2003.06.23
// Revised  N. Saoulidou 2004.03.05
// Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and 
// <a href="../CandTrackSR.html"> CandTrackSR Classes</a> (part of
// <a href="../index.html">The MINOS Class User Guide</a>)End_Html

#include <cassert>

#include "RawData/RawDigit.h"
#include "RawData/RawHeader.h"
#include "RawData/RawRecord.h"
#include "RawData/RawChannelId.h"
#include "RawData/RawDaqSnarlHeader.h"
#include "RawData/RawDaqHeaderBlock.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawVarcErrorInTfBlock.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "Algorithm/AlgConfig.h"
#include "Candidate/CandContext.h"
#include "CandData/CandHeader.h"
#include "CandTrackSR/AlgTrackSRList.h"
#include "CandTrackSR/CandTrackSR.h"
#include "CandTrackSR/CandTrackSR.h"
#include "CandTrackSR/CandTrackSRHandle.h"
#include "CandTrackSR/CandTrackSRList.h"
#include "CandTrackSR/CandTrackSRListHandle.h"
#include "CandTrackSR/HoughTrackSR.h"
#include "CandTrackSR/HoughViewSR.h"
#include "CandTrackSR/TrackClusterSR.h"
#include "CandTrackSR/Track2DSR.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStrip.h"
#include "Conventions/Mphysical.h"
#include "Conventions/Munits.h"
#include "Conventions/PlaneView.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Navigation/XxxItr.h"
#include "RecoBase/ArrivalTime.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/LinearFit.h"
#include "RecoBase/PropagationVelocity.h"
#include "Validity/VldContext.h"
#include "CandDigit/CandDeMuxDigitHandle.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "DataUtil/PlaneOutline.h"
#include "TMath.h"

//make a NavKey to select only U View strips that are not vetoed by demux
NavKey KeyOnUViewVetoCharge( const CandStripHandle *csh ){
  return (csh->GetPlaneView() == PlaneView::kU && !csh->GetDemuxVetoFlag());
}

//calorimeter ends with plane 120, so take all less than 121
NavKey KeyOnUViewVetoChargeNear( const CandStripHandle *csh ){
  return (csh->GetPlaneView() == PlaneView::kU && !csh->GetDemuxVetoFlag() && csh->GetPlane()<121);
}

//make a NavKey to select only V View strips that are not vetoed by demux
NavKey KeyOnVViewVetoCharge( const CandStripHandle *csh ){
  return (csh->GetPlaneView() == PlaneView::kV && !csh->GetDemuxVetoFlag());
}
NavKey KeyOnVViewVetoChargeNear( const CandStripHandle *csh ){
  return (csh->GetPlaneView() == PlaneView::kV && !csh->GetDemuxVetoFlag() && csh->GetPlane()<121);
}


//make a NavKey to sort Clusters by plane number
NavKey KeyOnClusterPlane( const TrackClusterSR *tc){
//   Int_t iplane = tc->GetPlane();
//   Float_t tpos = tc->GetMinTPos();
//   Int_t itpos = static_cast<Int_t>(tpos);
//   Float_t time = tc->GetBegTime();
//   Int_t itime = static_cast<Int_t>(time*.1/18.7);
  
//   Int_t navkey = (iplane<<22) | (itpos<<11) | itime;
  
//   return navkey;

  return tc->GetPlane();
}

//make a NavKey to select only U View Clusters
NavKey KeyOnUViewClusters( const TrackClusterSR *tc){
  return tc->GetPlaneView() == PlaneView::kU;
}

//make a NavKey to select only V View Clusters
NavKey KeyOnVViewClusters( const TrackClusterSR *tc){
  return tc->GetPlaneView() == PlaneView::kV;
}

//make a NavKey to select only U View tracks
NavKey KeyOnUViewTracks( const Track2DSR *tc){
  return tc->GetPlaneView() == PlaneView::kU;
}

//make a NavKey to select only V View tracks
NavKey KeyOnVViewTracks( const Track2DSR *tc){
  return tc->GetPlaneView() == PlaneView::kV;
}

//define the NavKey to select u view
static NavKey KeyOnUCluster( const TrackClusterSR *tc){
  return tc->GetPlaneView() == PlaneView::kU;
}

//define the NavKey to select v view
static NavKey KeyOnVCluster( const TrackClusterSR *tc){
  return tc->GetPlaneView() == PlaneView::kV;
}


//define the NavKey to select u view, non-wide clusters
static NavKey KeyOnUClusterNotWide( const TrackClusterSR *tc){
  return !tc->IsWide() && tc->GetPlaneView() == PlaneView::kU;
}

//define the NavKey to select v view, non-wide clusters
static NavKey KeyOnVClusterNotWide( const TrackClusterSR *tc){
  return !tc->IsWide() && tc->GetPlaneView() == PlaneView::kV;
}

//define the NavKey to select u view, non-wide clusters from 
//non-showerlike planes
static NavKey KeyOnUClusterNotWideSL( const TrackClusterSR *tc){
  return !tc->IsWide() && tc->GetPlaneView() == PlaneView::kU 
          && !tc->InShowerLikePlane();
}

//define the NavKey to select v view, non-wide clusters from 
//non-showerlike planes
static NavKey KeyOnVClusterNotWideSL( const TrackClusterSR *tc){
  return !tc->IsWide() && tc->GetPlaneView() == PlaneView::kV 
          && !tc->InShowerLikePlane();
}

//define the NavKey to select all u view clusters, from non-showerlike 
//planes
static NavKey KeyOnUClusterNotSL( const TrackClusterSR *tc){
  return tc->GetPlaneView() == PlaneView::kU && !tc->InShowerLikePlane();
}

//define the NavKey to select all v view clusters from non-showerlike 
//planes
static NavKey KeyOnVClusterNotSL( const TrackClusterSR *tc){
  return tc->GetPlaneView() == PlaneView::kV && !tc->InShowerLikePlane();
}


ClassImp(AlgTrackSRList)


//______________________________________________________________________
CVSID("$Id: AlgTrackSRList.cxx,v 1.109 2007/02/04 06:03:24 rhatcher Exp $");

//______________________________________________________________________
AlgTrackSRList::AlgTrackSRList()
{
  for(Int_t i = 0; i<1000; ++i){
    fMapIsWide[i] = 0;
  }
}

//______________________________________________________________________
AlgTrackSRList::~AlgTrackSRList()
{
}

//______________________________________________________________________
void AlgTrackSRList::RunAlg(AlgConfig &ac, CandHandle &ch,CandContext &cx)
{
  MSG("Alg", Msg::kDebug) << "Starting AlgTrackSRList::RunAlg()" << endl;

  assert(cx.GetDataIn());

  if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
    return;
  }

  CandTrackSRListHandle &cth = dynamic_cast<CandTrackSRListHandle &>(ch);

  const CandSliceListHandle *slicelist = 0;
  const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
  for (Int_t i=0; i<=cxin->GetLast(); i++) {
    TObject *tobj = cxin->At(i);
    if (tobj->InheritsFrom("CandSliceListHandle")) {
      slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
    }
  }
  fDetector = Detector::kUnknown;
  fSimFlag = SimFlag::kUnknown;

  if (!slicelist) {
    MSG("error",Msg::kError) <<
      "CandSliceListHandle missing\n";
  }
  else {
    fDetector = slicelist->GetVldContext()->GetDetector();
    fSimFlag = slicelist->GetVldContext()->GetSimFlag();
    MSG("TrackSR", Msg::kDebug) << "event number " 
                                 << slicelist->GetCandRecord()->GetCandHeader()->GetSnarl() 
                                 << endl;
  }


// Create Candcontext
  CandContext cxx(this,cx.GetMom());

// Get singleton instance of AlgFactory
  AlgFactory &af = AlgFactory::GetInstance();

  fParmMaxNStrip = ac.GetInt("MaxNStrip");
  fParmMaxNExtraStrip = ac.GetInt("MaxNExtraStrip");
  fParmTrk2DPlnEnd = ac.GetInt("Trk2DPlnEnd");
  fParmTrk2DWin0 = ac.GetDouble("Trk2DWin0");
  fParmHitNTime = ac.GetDouble("HitNTime");
  fParmHitTime = ac.GetDouble("HitTime");
  fParmTrk2DNSkip = ac.GetInt("Trk2DNSkip");
  fParmTrk2DNSkipHit = ac.GetInt("Trk2DNSkipHit");
  fParmTrk2DAlpha = ac.GetDouble("Trk2DAlpha");
  fParmTrk2DNSkipRemove = ac.GetInt("Trk2DNSkipRemove");
  fParmTrkNPlaneGoodDir = ac.GetInt("TrkNPlaneGoodDir");
  fParmTrk2DDirCosNSigma = ac.GetDouble("Trk2DDirCosNSig") ;
  fParmTrk2DDirCosError2 = ac.GetDouble("Trk2DDirCosError")*ac.GetDouble("Trk2DDirCosError");
  fParmTrk2DEndPlaneDiff = ac.GetInt("Trk2DEndPlaneDiff");
  fParmMinSingleHit = ac.GetInt("MinSingleHit");
  fParmSingleHitDef = ac.GetInt("SingleHitDef");
  fParmTrk2DNHough0 = ac.GetInt("Trk2DNHough0");
  fParmTrk2DLinA0 = ac.GetDouble("Trk2DLinA0");
  fParmTrk2DLinB0 = ac.GetDouble("Trk2DLinB0");
  fParmTrk2DNHough = ac.GetInt("Trk2DNHough");
  fParmTrk2DLinA = ac.GetDouble("Trk2DLinA");
  fParmTrk2DLinB = ac.GetDouble("Trk2DLinB");
  fParmTrk2DNSigmaA = ac.GetDouble("Trk2DNSigmaA");
  fParmTrk2DNSeed = ac.GetInt("Trk2DNSeed");
  fParmTrk2DMaxResid = ac.GetDouble("Trk2DMaxResid");
  fParmTrk2DNContiguous = ac.GetInt("Trk2DNContiguous");
  fParmTrk2DHitFraction = ac.GetDouble("Trk2DHitFraction");
  fParmTrk2DTwinMatchFraction = ac.GetDouble("Trk2DTwinMatchFraction");
  fParmTrk2DSubsetNHit = ac.GetInt("Trk2DSubsetNHit");
  fParmTrk2DSubsetDHit = ac.GetInt("Trk2DSubsetDHit");
  fParmDiffViewBegPlaneMatch = ac.GetInt("DiffViewBegPlaneMatch");
  fParmDiffViewEndPlaneMatch = ac.GetInt("DiffViewEndPlaneMatch");
  fParmDiffViewPlaneOverlap = ac.GetDouble("DiffViewPlaneOverlap");
  fParmDiffViewTimeMatch = ac.GetDouble("DiffViewTimeMatch");
  fParmTrk3DTwinFrac = ac.GetDouble("Trk3DTwinFrac");
  fParmTrkClsNSkip = ac.GetInt("TrkClsNSkip");
  fParmTrk3DTrackMax = ac.GetInt("Trk3DTrackMax");
  fParmIsCosmic = ac.GetInt("IsCosmic");
  fParmMinStripPulseHeight = ac.GetDouble("MinStripPulseHeight");
  fParmMinClusterPulseHeight = ac.GetDouble("MinClusterPulseHeight");
  fParmUseWideClusters = ac.GetInt("UseWideClusters");
  fParmUseShowerLikePlanes = ac.GetInt("UseShowerLikePlanes");
  fParmMaxTimingResid = ac.GetInt("MaxTimingResid")*Munits::ns;
  fParmMisalignmentError = ac.GetInt("MisalignmentError")*Munits::mm;
  fParmExtrapError = ac.GetDouble("ExtrapError");

  fParmHoughMinBin =  fParmTrk2DNHough0 ;
  fParmHoughMinFrac = 0.5;
  fParmHoughMinFracZoom = .75;
  fParmHoughMinInterBinSize = .02;

  int EnableSpectFind = 1;
  EnableSpectFind = ac.GetInt("EnableSpectFind");

  const CandRecord *candrec = cx.GetCandRecord();
  assert(candrec);
  const VldContext *vldcptr = candrec->GetVldContext();
  assert(vldcptr);
  fVldc = *vldcptr; 
  
  UgliGeomHandle ugh(fVldc);

 // Get an AlgHandle to AlgTrackSR with proper AlgConfig
  const char *pTrackAlgConfig = 0;
  ac.Get("TrackAlgConfig",pTrackAlgConfig);
  AlgHandle ah = af.GetAlgHandle("AlgTrackSR",pTrackAlgConfig);


  // set new strip and digit list handles to null at this point
  // if we have near detector spectrometer tracking these will get set 
  // during process of the first slice
  fCDLHnew=0;
  fCSLHnew=0;
  fCSlcLHnew=0;

  //get an iterator over the slice list passed into the algorithm
  //and loop over all slices
  Int_t islice=0;

  CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
  while (CandSliceHandle *slice = sliceItr()) {
    MSG("TrackSR", Msg::kDebug) << " ******* Slice ********* " << islice <<  endl;
  
    islice++;

    //make variables to keep track of the event extent in 
    //the different views.  also keep track of the first and 
    //last planes hit in the event and the number of planes 
    //hit in a given view.
    fSliceVtxUPlane = 1000;
    fSliceVtxVPlane = 1000;
    fSliceEndUPlane = -1;
    fSliceEndVPlane = -1;
    fSliceVtxPlane = 1000;
    fSliceEndPlane = -1;
    fSliceNUPlanes = 0;
    fSliceNVPlanes = 0;
    fSliceNPlanes = 0;
    
    //zero the map of wide planes for this slice
    for(Int_t i = 0; i<1000; ++i){
      fMapIsWide[i] = 0;
    }
    
    //get a couple of iterators over the strips in the slice    
    CandStripHandleItr uStripIter(slice->GetDaughterIterator());        
    CandStripHandleItr vStripIter(slice->GetDaughterIterator());        
    CandStripHandleItr allStripIter(slice->GetDaughterIterator());      
        
     //sort the strips using the 2 function sort in navigation  
        uStripIter.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime);
        vStripIter.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime);
        allStripIter.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime);

    /* for now leave in the old sort just in case       
    //set the key functions to sort the strips by plane, strip and time
    CandStripHandleKeyFunc *uStripKF = uStripIter.CreateKeyFunc();
    uStripKF->SetFun(CandStripHandle::KeyFromPlaneStripTime);
    uStripIter.GetSet()->AdoptSortKeyFunc(uStripKF);
    uStripKF = 0;

    CandStripHandleKeyFunc *vStripKF = vStripIter.CreateKeyFunc();
    vStripKF->SetFun(CandStripHandle::KeyFromPlaneStripTime); 
    vStripIter.GetSet()->AdoptSortKeyFunc(vStripKF);
    vStripKF = 0;

    CandStripHandleKeyFunc *allStripKF = allStripIter.CreateKeyFunc();
    allStripKF->SetFun(CandStripHandle::KeyFromPlaneStripTime);
    allStripIter.GetSet()->AdoptSortKeyFunc(allStripKF);
    allStripKF = 0;
    */
    allStripIter.ResetFirst();
   
    // a loop to check that the sorting worked
    /*
    CandStripHandle *strip = 0;
    while( (strip = allStripIter()) ){
      MSG("TrackSR", Msg::kDebug) << "plane " << strip->GetPlane() << " strip " 
                                  << strip->GetStrip() << " " 
                                  << (Int_t)strip->GetPlaneView() << " " 
                                  << (Int_t)PlaneView::kX << " " 
                                  << (Int_t)PlaneView::kY << " " 
                                  << (Int_t)PlaneView::kU << " " 
                                  << (Int_t)PlaneView::kV << " " 
                                  << endl;
    } 
    */

    //now adopt a select key function to only get those strips
    //in each view that were not vetoed 
    CandStripHandleKeyFunc *uSelectKF = uStripIter.CreateKeyFunc();
    if(fDetector!=Detector::kNear){
      uSelectKF->SetFun(KeyOnUViewVetoCharge);
    }
    else{
      uSelectKF->SetFun(KeyOnUViewVetoChargeNear);
    }

    uStripIter.GetSet()->AdoptSelectKeyFunc(uSelectKF);
    uSelectKF = 0;

    CandStripHandleKeyFunc *vSelectKF = vStripIter.CreateKeyFunc();
    if(fDetector!=Detector::kNear){
      vSelectKF->SetFun(KeyOnVViewVetoCharge);
    }
    else{
      vSelectKF->SetFun(KeyOnVViewVetoChargeNear);
    }
    vStripIter.GetSet()->AdoptSelectKeyFunc(vSelectKF);
    vSelectKF = 0;      
    
    uStripIter.ResetFirst();
    vStripIter.ResetFirst();
    allStripIter.ResetFirst();

    //find the vtx and end planes, as well
    //as the total number of planes in each view
    Find2DTrackEndPoints(uStripIter, 
                         fSliceVtxUPlane, 
                         fSliceEndUPlane, 
                         fSliceNUPlanes);
    Find2DTrackEndPoints(vStripIter, 
                         fSliceVtxVPlane, 
                         fSliceEndVPlane, 
                         fSliceNVPlanes);
    MSG("TrackSR", Msg::kDebug) << "planes in slice " << fSliceNVPlanes << " " 
                                << fSliceNUPlanes << " " << fParmTrk2DNSeed << endl;

    if(fSliceNUPlanes<fParmTrk2DNSeed || fSliceNVPlanes<fParmTrk2DNSeed){
      MSG("TrackSR", Msg::kDebug) << "too few planes in slice" << endl;
      continue;
    }

    fSliceNPlanes = fSliceNUPlanes + fSliceNVPlanes;

    //get the rough tracking information from the HoughViewSR objects
    //these give you an idea of the slope and intercept for the 2D tracks
    //in the event, as well as the number of tracks to expect in the event
    //the SetStripList method returns true if there are enough planes in
    //each view that pass the requirements for estimating using the hough transform

    HoughViewSR houghUView;
    houghUView.SetPeakMin(fParmHoughMinBin);
    houghUView.SetPeakMinFrac(fParmHoughMinFrac);
    houghUView.SetPeakMinFracZoom(fParmHoughMinFracZoom);
    houghUView.SetMinInterBinSize(fParmHoughMinInterBinSize);

    HoughViewSR houghVView;
    houghVView.SetPeakMin(fParmHoughMinBin);
    houghVView.SetPeakMinFrac(fParmHoughMinFrac);
    houghVView.SetPeakMinFracZoom(fParmHoughMinFracZoom);
    houghVView.SetMinInterBinSize(fParmHoughMinInterBinSize);

    //quit trying to track if you dont have enough 

    Int_t maxUTrack = 0;
    Int_t maxVTrack = 0;
    
    Double_t uTrackSlope = 0.;
    Double_t uTrackIntercept = 0.;
    Double_t vTrackSlope = 0.;
    Double_t vTrackIntercept = 0.;

    //get the slopes and intercepts for the different views
 
    TObjArray *utrackClusterList = new TObjArray(1,0);   
    //fill the cluster list one view at a time

    MakeTrackClusters(utrackClusterList, uStripIter, cth, 1, 1);
    TrackClusterSR *tc = 0;
    Double_t xfit[500];
    Double_t yfit[500];
    Double_t wfit[500];
    Double_t parm[2];
    Double_t eparm[2] = {0.,0.};
    Double_t wtmp;
    Int_t ifit=0;
    for(Int_t i=0; i<=utrackClusterList->GetLast(); ++i){
      tc = dynamic_cast<TrackClusterSR*>(utrackClusterList->At(i));
      if(ifit<500){
        Bool_t use=true;
        for(Int_t j=0; j<=utrackClusterList->GetLast(); ++j){
          TrackClusterSR * tc2 = dynamic_cast<TrackClusterSR*>(utrackClusterList->At(j));
          if(i!=j  && tc->GetPlane()==tc2->GetPlane()){
            use=false;
            break;
          }
        }
        if(use){
          xfit[ifit] = tc->GetZPos();
          yfit[ifit] = tc->GetTPos();
          wtmp = 0.288*(tc->GetMaxTPos()-tc->GetMinTPos());
          wfit[ifit]=1.0/(wtmp*wtmp);
          ifit++;
        }
      }
    }
    if(ifit>2){
      LinearFit::Weighted(ifit,xfit,yfit,wfit,parm,eparm);
      uTrackSlope=parm[1];
      uTrackIntercept=parm[0];
    }

    houghUView.SetClusterList(utrackClusterList, fParmMinStripPulseHeight); 

    houghUView.Iterate();
    MSG("TrackSR", Msg::kDebug) << "u tracks " << houghUView.GetNTrack() << endl;
    if( houghUView.GetNTrack() > 0 ){
      maxUTrack = houghUView.GetNTrack();
      const HoughTrackSR *ht = houghUView.GetHoughTrack(houghUView.GetLargestTrackIndex());
      uTrackSlope = ht->GetPeakSlope(ht->GetLargestPeakIndex());
      uTrackIntercept = ht->GetPeakIntercept(ht->GetLargestPeakIndex());
    }

    TObjArray *vtrackClusterList = new TObjArray(1,0);   
    //fill the cluster list one view at a time
    MakeTrackClusters(vtrackClusterList, vStripIter, cth, 1, 1);
   
    ifit=0;
    for(Int_t i=0; i<=vtrackClusterList->GetLast(); ++i){
      tc = dynamic_cast<TrackClusterSR*>(vtrackClusterList->At(i));
      if(ifit<500){
        Bool_t use=true;
        for(Int_t j=0; j<=vtrackClusterList->GetLast(); ++j){
          TrackClusterSR * tc2 = dynamic_cast<TrackClusterSR*>(vtrackClusterList->At(j));
          if(i!=j  && tc->GetPlane()==tc2->GetPlane()){
            use=false;
            break;
          }
        }
        if(use){
          xfit[ifit] = tc->GetZPos();
          yfit[ifit] = tc->GetTPos();
          wtmp = 0.288*(tc->GetMaxTPos()-tc->GetMinTPos());
          wfit[ifit]=1.0/(wtmp*wtmp);
          ifit++;
        }
      }
    }
    if(ifit>2){
      LinearFit::Weighted(ifit,xfit,yfit,wfit,parm,eparm);
      vTrackSlope=parm[1];
      vTrackIntercept=parm[0];
    }
    houghVView.SetClusterList(vtrackClusterList, fParmMinStripPulseHeight); 
    houghVView.Iterate();
    MSG("TrackSR", Msg::kDebug) << "v tracks " << houghVView.GetNTrack() << endl;
    if( houghVView.GetNTrack() > 0 ){
      maxVTrack = houghVView.GetNTrack();
      const HoughTrackSR *ht = houghVView.GetHoughTrack(houghVView.GetLargestTrackIndex());
      vTrackSlope = ht->GetPeakSlope(ht->GetLargestPeakIndex());
      vTrackIntercept = ht->GetPeakIntercept(ht->GetLargestPeakIndex());
    }


 // delete track clusters
    for(Int_t i=0; i<=utrackClusterList->GetLast(); ++i){
      tc = dynamic_cast<TrackClusterSR*>(utrackClusterList->At(i));
      delete tc;
    }
    delete utrackClusterList;

    for(Int_t i=0; i<=vtrackClusterList->GetLast(); ++i){
      tc = dynamic_cast<TrackClusterSR*>(vtrackClusterList->At(i));
      delete tc;
    }
    delete vtrackClusterList;
    
    
    //get the slope of the track wrt the z direction, dz/ds
    fSliceDzDs = 1./TMath::Sqrt(1.+uTrackSlope*uTrackSlope+vTrackSlope*vTrackSlope);

    MSG("TrackSR",Msg::kDebug) << "Hough: # tracks = " << maxUTrack 
                                 << " " << maxVTrack 
                                 << endl
                                 << "Hough: slope = " << uTrackSlope 
                                 << " " << vTrackSlope 
                                 << endl
                                 << "Hough: intercept = " << uTrackIntercept << " " 
                                 << vTrackIntercept << endl;
    
    if(uTrackSlope==-1)uTrackSlope=0;
    if(vTrackSlope==-1)vTrackSlope=0;

    fHoughUSlope = uTrackSlope;
    fHoughVSlope = vTrackSlope;
    fHoughUIntercept = uTrackIntercept;
    fHoughVIntercept = vTrackIntercept;


    //set the initial fit direction with respect to z 
    //1 = forward, -1 = backwards
    Int_t idir = 1; 
    //if this is a beam file, just set the direction to -1 because
    //it is easier to work backwards to the vertex
    if(!fParmIsCosmic) idir = -1; 

    //figure out where the vtx and end planes for the event
    fSliceVtxPlane = TMath::Min(fSliceVtxUPlane,fSliceVtxVPlane);
    fSliceEndPlane = TMath::Max(fSliceEndUPlane,fSliceEndVPlane);

    //determine whether it is worth it to do the timing for the event
    //make sure at least 1 track from the hough transform in each view
    if(maxUTrack>0 && maxVTrack>0){
      allStripIter.ResetFirst();
      //check and see if you have the direction right using the
      //timing in the event
      if(FindTimingDirection(allStripIter, uTrackSlope, uTrackIntercept,
                             vTrackSlope, vTrackIntercept)<0 ){
        if(fParmIsCosmic){
          idir = -1;
          fSliceEndPlane = TMath::Min(fSliceVtxUPlane,fSliceVtxVPlane);
          fSliceVtxPlane = TMath::Max(fSliceEndUPlane,fSliceEndVPlane);
        }
      }
    }//end if enough tracks to do the timing
    
    allStripIter.ResetFirst();

    Int_t uExpectedStrips = fParmMaxNStrip-fParmMaxNExtraStrip;
    Int_t vExpectedStrips = fParmMaxNStrip-fParmMaxNExtraStrip;

    MSG("TrackSR", Msg::kDebug) << "max u expected strips " << uExpectedStrips
                               << " max v expected strips " << vExpectedStrips 
                               << endl;

    //the number of expected strips is based on the slope of the track in
    //tpos vs z /4.108 where 4.108 is the width of a strip in tpos
    //so the value is how many strips you think will be hit/plane
    Int_t expectedStrips = (Int_t)(TMath::Abs(uTrackSlope)/4.108)+1;
    
    //only change the # of expected strips if you did the hough tracking
    //and it returned with at least 1 track in each view
    if(expectedStrips<uExpectedStrips && maxUTrack>0) 
      uExpectedStrips = expectedStrips;
    
    expectedStrips = (Int_t)(TMath::Abs(vTrackSlope)/4.108)+1;
    
    if(expectedStrips<vExpectedStrips && maxVTrack>0) 
      vExpectedStrips = expectedStrips;

    MSG("TrackSR", Msg::kDebug) << "u expected strips " << uExpectedStrips
                               << " v expected strips " << vExpectedStrips 
                               << " u track slope " << uTrackSlope 
                               << " v track slope " << vTrackSlope << endl;
      
    //now we can make track clusters.  declare an array to hold the cluster
    //objects, but fill it in the MakeTrackClusters() method
    TObjArray *trackClusterList = new TObjArray(1,0);
   
    //fill the cluster list one view at a time
    MakeTrackClusters(trackClusterList, uStripIter, cth, uExpectedStrips, idir);
    MakeTrackClusters(trackClusterList, vStripIter, cth, vExpectedStrips, idir);
      
    //an iterator over all clusters in the event
    TrackClusterSRItr masterClusterItr(trackClusterList);
    TrackClusterSRKeyFunc *masterClusterKf = masterClusterItr.CreateKeyFunc();
    masterClusterKf->SetFun(KeyOnClusterPlane);
    masterClusterItr.GetSet()->AdoptSortKeyFunc(masterClusterKf);
    masterClusterKf = 0;
    
    //clusterItr has all clusters in a view that pass the given cuts
    TrackClusterSRItr clusterItr(trackClusterList);
    TrackClusterSRKeyFunc *clusterKf = clusterItr.CreateKeyFunc();
    clusterKf->SetFun(KeyOnClusterPlane);
    clusterItr.GetSet()->AdoptSortKeyFunc(clusterKf);
    clusterKf = 0;
 
 
    //planeClusterItr will keep track of the clusters in a given plane
    TrackClusterSRItr planeClusterItr(trackClusterList);
    TrackClusterSRKeyFunc *planeClusterKf = planeClusterItr.CreateKeyFunc();
    planeClusterKf->SetFun(KeyOnClusterPlane);
    planeClusterItr.GetSet()->AdoptSortKeyFunc(planeClusterKf);
    planeClusterKf = 0;
    
    //fill array of ints to ints where the index is the 
    //plane number and the value is 0 for planes
    //containing only good clusters, and 1 for planes with
    //wide clusters
    
    clusterItr.ResetFirst();
    TrackClusterSR *cluster = 0;
    while( (cluster = clusterItr()) ){
      if( cluster->IsWide() ) fMapIsWide[cluster->GetPlane()] = (Int_t)cluster->IsWide();

      MSG("TrackSR", Msg::kDebug) << "cluster on plane " << cluster->GetPlane() 
                                 << " tpos " << cluster->GetTPos() << " fMapIsWide "
                                 << fMapIsWide[cluster->GetPlane()] << endl;
    }
    clusterItr.ResetFirst();
    
    //check to see how many planes you might expect in a track
    //by looking for gaps in the planes hit.
 
    Int_t nmax3dtrk = 0;
    int trk2dmin = fParmTrk2DNSeed;
   
    
    //
    // 2D tracking
    //
    
    
    //loop over each view separately to do the 2d tracking
    Int_t viewCtr = 0;
    Double_t houghSlope = 0.;
    Double_t houghIntercept = 0.;
    Int_t nmaxtrack = 0;
    Int_t sliceNPlane = 0;
    Int_t sliceVtxPlane = 0;
    Int_t sliceEndPlane = 0;
    
    //create a TObjArray to hold the 2d tracks we find
    TObjArray *trackList = new TObjArray(1,0);
    TObjArray *seedClusters = new TObjArray(1,0);
    
    while(viewCtr<2){
      
      //first we want to select only those clusters for this view
      //we dont want to use any wide clusters when doing the track finding
      //depending on whether we will allow clusters from showerlike planes
      //we either use KeyOnXClusterNotWide or KeyOnXClusterNotWideSL
      //where the X is either U or V and NotWide excludes only wide
      //clusters, NotWideSL excludes wide clusters and those from showerlike
      //planes.  if we allow wide clusters, then use KeyOnXCluster to just
      //select a view
      
      //the master cluster list has all clusters that pass the min 
      //pulseheight cut.  if wide clusters are allowed it will have 
      //those too.  it will always have clusters from showerlike planes.
      //so it either uses KeyOnU(V)Cluster or KeyOnU(V)ClusterNotWide
      
      TrackClusterSRKeyFunc *selectKF = clusterItr.CreateKeyFunc();

      TrackClusterSRKeyFunc *planeSelectKF = planeClusterItr.CreateKeyFunc();
      TrackClusterSRKeyFunc *masterSelectKF = masterClusterItr.CreateKeyFunc();

      //u view first
      if(viewCtr == 0){
        houghSlope = uTrackSlope;
        houghIntercept = uTrackIntercept;
        nmaxtrack = maxUTrack;
        sliceNPlane = fSliceNUPlanes;
        sliceVtxPlane = fSliceVtxUPlane;
        sliceEndPlane = fSliceEndUPlane;
        
        if(fParmUseWideClusters) masterSelectKF->SetFun(KeyOnUCluster);
        else masterSelectKF->SetFun(KeyOnUClusterNotWide);
        
        if(fParmUseWideClusters && fParmUseShowerLikePlanes){
          selectKF->SetFun(KeyOnUCluster);
          planeSelectKF->SetFun(KeyOnUCluster);
        }
        else if(fParmUseWideClusters && !fParmUseShowerLikePlanes){
          selectKF->SetFun(KeyOnUClusterNotSL);
          planeSelectKF->SetFun(KeyOnUClusterNotSL);
        }
        else if(!fParmUseWideClusters && fParmUseShowerLikePlanes){
          selectKF->SetFun(KeyOnUClusterNotWide);
          planeSelectKF->SetFun(KeyOnUClusterNotWide);
        }
        else if(!fParmUseWideClusters && !fParmUseShowerLikePlanes){
          selectKF->SetFun(KeyOnUClusterNotWideSL);
          planeSelectKF->SetFun(KeyOnUClusterNotWideSL);
        }
      }//end if 0 view
      else if( viewCtr == 1){
        houghSlope = vTrackSlope;
        houghIntercept = vTrackIntercept;
        nmaxtrack = maxVTrack;
        sliceNPlane = fSliceNVPlanes;
        sliceVtxPlane = fSliceVtxVPlane;
        sliceEndPlane = fSliceEndVPlane;
        
        if(fParmUseWideClusters) masterSelectKF->SetFun(KeyOnVCluster);
        else masterSelectKF->SetFun(KeyOnVClusterNotWide);
        
        if(fParmUseWideClusters && fParmUseShowerLikePlanes){
          selectKF->SetFun(KeyOnVCluster);
          planeSelectKF->SetFun(KeyOnVCluster);
        }
        else if(fParmUseWideClusters && !fParmUseShowerLikePlanes){
          selectKF->SetFun(KeyOnVClusterNotSL);
          planeSelectKF->SetFun(KeyOnVClusterNotSL);
        }
        else if(!fParmUseWideClusters && fParmUseShowerLikePlanes){
          selectKF->SetFun(KeyOnVClusterNotWide);
          planeSelectKF->SetFun(KeyOnVClusterNotWide);
        }
        else if(!fParmUseWideClusters && !fParmUseShowerLikePlanes){
          selectKF->SetFun(KeyOnVClusterNotWideSL);
            planeSelectKF->SetFun(KeyOnVClusterNotWideSL);
        }
      }//end if v view
      
      clusterItr.GetSet()->AdoptSelectKeyFunc(selectKF);
      planeClusterItr.GetSet()->AdoptSelectKeyFunc(planeSelectKF);
      masterClusterItr.GetSet()->AdoptSelectKeyFunc(masterSelectKF);
      selectKF = 0;
      planeSelectKF = 0;
      masterSelectKF = 0;
      
      TrackClusterSR *cluster;

      masterClusterItr.ResetFirst();
      
      //set the iterator direction to fill the vector of hit planes
      if(idir == 1){
        clusterItr.ResetFirst();
        planeClusterItr.ResetFirst();
      }
      else if(idir == -1){
        clusterItr.ResetLast();
        planeClusterItr.ResetLast();
      }
      
      //create a vector of ints to hold the numbers of 
      //the planes with acceptable clusters in them
      //only holds the numbers for  planes in the current view
      vector<Int_t> hitPlanes;
      Int_t prevPlane = 0;
      while( (cluster = clusterItr()) ){
        
        MSG("TrackSR", Msg::kDebug) << "cluster on plane " 
                                   << cluster->GetPlane() << " time = " 
                                   << cluster->GetBegTime() << endl;
        if(cluster->GetPlane() != prevPlane){
          hitPlanes.push_back(cluster->GetPlane());
          prevPlane = cluster->GetPlane();
        }
        
      }//end loop over clusters in the iterator
      
      FindTimingDirection(clusterItr);
      
      Int_t iterate = 0;
      bool isDone = false;
      Int_t clustersLeft = clusterItr.SizeSelect();
      Int_t begPlane = 0;
      Int_t loopTracks = 1;
      
      //do the following loop until the track is done or
      //there are no more valid clusters left
      //     while( !isDone && clustersLeft>0 && loopTracks>0){
      while( !isDone && clustersLeft>0){
        
        //set the iterator direction to fill the vector of hit planes
        if(idir == 1){
          clusterItr.ResetFirst();
          planeClusterItr.ResetFirst();
        }
        else if(idir == -1){
          clusterItr.ResetLast();
          planeClusterItr.ResetLast();
        }
        Track2DSR *track = 0;
        Int_t trackCtr = 0;     
        //reset the number of clusters left in this view
        clustersLeft = 0;
        
        //identify the seed clusters and put them into trackList
        loopTracks = FillTrackList(clusterItr, trackList, seedClusters,
                                   houghSlope, houghIntercept, idir, iterate, 
                                   begPlane);
        
        //get an iterator over the tracks
        Track2DSRItr trackItr(trackList);
        Track2DSRItr trackItr2(trackList);
        
        //program the selection key functions to get only tracks
        //in the current view
        Track2DSRKeyFunc *select1KF = trackItr.CreateKeyFunc();
        Track2DSRKeyFunc *select2KF = trackItr2.CreateKeyFunc();
        
        if(viewCtr == 0){
          select1KF->SetFun(KeyOnUViewTracks);
          select2KF->SetFun(KeyOnUViewTracks);
        }
        else if(viewCtr == 1){
          select1KF->SetFun(KeyOnVViewTracks);
          select2KF->SetFun(KeyOnVViewTracks);
        }
        trackItr.GetSet()->AdoptSelectKeyFunc(select1KF);
        trackItr2.GetSet()->AdoptSelectKeyFunc(select2KF);
        
        select1KF = 0;
        select2KF = 0;


        MSG("TrackSR",Msg::kDebug) << "after FillTrackList, total # tracks = " 
                                  << trackList->GetLast()+1 << "\n";

        //loop over the hit planes and tracks made in this iteration
        //to add clusters to the tracks
        AddClustersToTracks(trackItr, planeClusterItr, idir, iterate, 
                            hitPlanes,trk2dmin);


        MSG("TrackSR",Msg::kDebug) << " after add clusters, total # tracks = " 
                                  << trackList->GetLast()+1 << "\n";

        IdentifyBadTracks(trackItr, iterate, trk2dmin);

          trackItr.ResetFirst();
          trackCtr = 0;
          while( (track = trackItr()) ){

            MSG("TrackSR", Msg::kDebug) << "track " << trackCtr 
                                       << " is bad " << (Int_t)track->IsBad() << endl; 
            for(Int_t i = 0; i<=track->GetLast(); ++i){
              cluster = track->GetCluster(i);
              MSG("TrackSR", Msg::kDebug) << "\tTC plane = " << cluster->GetPlane()
                                         << " tpos = " << cluster->GetTPos() 
                                         << endl;
            }
            ++trackCtr;
          }     
        MarkUsedClusters(trackItr, clusterItr, iterate, nmaxtrack, houghSlope, 
                         houghIntercept);

 
        //subtract the # of tracks from the count of that are formed and then
        //removed in this loop
        loopTracks -= RemoveTrackSubsets(trackItr, trackItr2, trackList);
         
        MSG("TrackSR",Msg::kDebug) << "after remove subsets, total # tracks = " 
                                  << trackList->GetLast()+1 << "\n";
          
        //need to count the number of valid clusters left in the clusterItr
        clusterItr.ResetFirst();
        while( (cluster = clusterItr()) ){
          if(cluster->IsValid()){
            ++clustersLeft;
            MSG("TrackSR", Msg::kDebug) << "valid cluster on plane " 
                                       << cluster->GetPlane() << endl;
          }
        }
        
        //see if the tracks span the event
        // if cosmic mode and found tracks span event stop
        if(fParmIsCosmic && sliceVtxPlane>0 
           && sliceEndPlane>0 ){
          Bool_t isDoneU=false;
          Bool_t isDoneV=false;
          for(Int_t itrk=0; itrk<=trackList->GetLast(); itrk++){
            Track2DSR *track = dynamic_cast<Track2DSR*>(trackList->At(itrk));
            
            //if the number of clusters in the track is the same as 
            //the number of hit planes in the view, stop tracking
            if(track->GetLast()+1==sliceNPlane){
              if(track->GetPlaneView()==PlaneView::kU) isDoneU = true;
              if(track->GetPlaneView()==PlaneView::kV) isDoneV = true;
              isDone = isDoneU & isDoneV;
              if(isDone)break;
            }
          }//end loop over tracks
        }//end if a cosmic event 

        //another iteration down, increment the counter
        ++iterate;
      }//end loop to find all the tracks in the view
            

     //now fill in the gaps in the tracks for this view
      Track2DSRItr trackItr(trackList);
      FillInGaps(trackItr, masterClusterItr, idir);
      
      //clear the iterator select functions for the next view
      masterClusterItr.GetSet()->AdoptSelectKeyFunc(0);
      planeClusterItr.GetSet()->AdoptSelectKeyFunc(0);
      clusterItr.GetSet()->AdoptSelectKeyFunc(0);
      
      //increment the view cnt to do the next view
      ++viewCtr;
    }//end loop over each view for 2d tracking
    
    // merge tracks which are non-overlapping, and which are consistent
    MergeTracks(trackList,idir);
    
    //get iterators over the u and v view tracks
    Track2DSRItr uTrackItr(trackList);
    Track2DSRItr vTrackItr(trackList);
    
    MSG("TrackSR", Msg::kDebug) << "number of 2d tracks = " << trackList->GetLast()+1  << endl;
    
    CandSliceHandle * trackSlice=slice;
    
    // ********************************************
    // for near detector, add spectrometer hits to tracks which end near
    // plane 121. 
    // *******************************************
    if(EnableSpectFind){
      if(fDetector==Detector::kNear){
        CandSliceHandle *  dupSlice=slice->DupHandle();
        SpectrometerTracking(trackList, cth, slice, dupSlice, cx);
        // if spectrometer tracking has cloned slice, associate that slice
        // with track
        if(dupSlice->GetUidInt()!=slice->GetUidInt()){
          trackSlice = dynamic_cast<CandSliceHandle *>
            (fCSlcLHnew->FindDaughter(dupSlice));
        }
        delete dupSlice;
      }  
    }
    //make the 3d tracks out of the 2d tracks
    TObjArray *candTracks = new TObjArray(1,0);
 
    //program the selection key functions
    Track2DSRKeyFunc *selectUKF = uTrackItr.CreateKeyFunc();
    Track2DSRKeyFunc *selectVKF = vTrackItr.CreateKeyFunc();

    
    selectUKF->SetFun(KeyOnUViewTracks);
    selectVKF->SetFun(KeyOnVViewTracks);
    
    uTrackItr.GetSet()->AdoptSelectKeyFunc(selectUKF);
    vTrackItr.GetSet()->AdoptSelectKeyFunc(selectVKF);
    
    selectUKF = 0;
    selectVKF = 0;
  
    FormCandTrackSR(uTrackItr, vTrackItr, candTracks, ch, trackSlice, ah, cxx, idir);
    
    //get iterators over the cand tracks
    CandTrackSRHandleItr candTracks1(candTracks);
    CandTrackSRHandleItr candTracks2(candTracks);
    
    //get rid of possible twin tracks
    DeleteTwinTracks(candTracks1, candTracks2, candTracks,
                     maxUTrack, maxVTrack,nmax3dtrk, ch);
    

    //free up the memory used by trackList and trackClusterList
    CleanUp(trackList, trackClusterList);    
    delete trackClusterList;
    delete trackList;
    delete seedClusters;
    delete candTracks;
    
  }//end loop over slices passed in

  // remove strips from striplist in spectrometer
  // if not on any track

  if(EnableSpectFind){
    if(fDetector==Detector::kNear){
      RemoveUnusedSpectStrips(cth,fCSLHnew);
      RemoveStripsInSlice(cth,fCSlcLHnew);
    }
  }

  if(fCDLHnew){
    CandRecord *crec = cx.GetCandRecord();
    assert(crec);
    fCDLHnew->SetName("canddigitlist");
    crec->SecureCandHandle(*fCDLHnew);
    delete fCDLHnew;
  }
  if(fCSlcLHnew){
    CandRecord *crec = cx.GetCandRecord();
    assert(crec);
    fCSlcLHnew->SetName("CandSliceList");
    crec->SecureCandHandle(*fCSlcLHnew);
    delete fCSlcLHnew;
  }
  if(fCSLHnew){ 
    CandRecord *crec = cx.GetCandRecord();
    assert(crec);
    fCSLHnew->SetName("CandStripList");
    crec->SecureCandHandle(*fCSLHnew);
    delete fCSLHnew;
  }

  MSG("TrackSR", Msg::kDebug) << "done with AlgTrackSRList" << endl;

  return;
}//end Run Alg


//-------------------------------------------
void AlgTrackSRList::RemoveUnusedSpectStrips(CandTrackSRListHandle & cth, 
                                            CandStripListHandle *striplist)
{
  if(!fDetector==Detector::kNear)return;
  if(striplist){
    CandStripHandleItr stripItr(striplist->GetDaughterIterator());
    for (CandStripHandle *strip = stripItr(); strip ; strip = stripItr()) {
      if(strip->GetPlane()>120){
        CandDigitHandleItr digitItr1(strip->GetDaughterIterator());
        CandDigitHandle *digit1 = dynamic_cast<CandDigitHandle*>(digitItr1());
        Bool_t stripIsOnTrack=false;
        Bool_t samePixel=false;
        TObjArray *tclist = cth.GetTrackClusterList();
        for (int i=0; i<=tclist->GetLast(); i++) {
          TrackClusterSR *tc = dynamic_cast<TrackClusterSR*>(tclist->At(i));
          if(tc->GetPlane()>120 ){
            const TObjArray *trackstriplist = tc->GetStripList();
            for (int j=0; j<=trackstriplist->GetLast(); j++) {
              CandStripHandle * trackstrip = dynamic_cast<CandStripHandle*>(trackstriplist->At(j));
              CandDigitHandleItr digitItr2(trackstrip->GetDaughterIterator());
              CandDigitHandle *digit2 = dynamic_cast<CandDigitHandle*>(digitItr2());

              if(trackstrip->IsEqual(strip)){
                stripIsOnTrack=true;
              }
              if(digit1->GetChannelId()==digit2->GetChannelId()){
                samePixel=true;
              }
            }
          }
        }
        if(!stripIsOnTrack && samePixel){
          striplist->RemoveDaughter(strip);
          
        }
      }
    }
  }
}

//-------------------------------------------
void AlgTrackSRList::RemoveStripsInSlice(CandTrackSRListHandle &cth,
                                         CandSliceListHandle * slicelist)
{
 
  if(!fDetector==Detector::kNear)return;
  if(slicelist){
    CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
    for (CandSliceHandle *slice = sliceItr(); slice ; slice = sliceItr()) {
      CandSliceHandle *  dupSlice=slice->DupHandle(); 
      CandStripHandleItr slicestripItr(dupSlice->GetDaughterIterator());
      for (CandStripHandle *strip = slicestripItr(); strip ; strip = slicestripItr()) {
        if(strip->GetPlane()>120){
          CandDigitHandleItr digitItr1(strip->GetDaughterIterator());
          CandDigitHandle *digit1 = dynamic_cast<CandDigitHandle*>(digitItr1());
          Bool_t stripIsOnTrack=false;
          Bool_t samePixel=false;  

        TObjArray *tclist = cth.GetTrackClusterList();
        for (int i=0; i<=tclist->GetLast(); i++) {
          TrackClusterSR *tc = dynamic_cast<TrackClusterSR*>(tclist->At(i));
          if(tc->GetPlane()>120 ){
            const TObjArray *trackstriplist = tc->GetStripList();
            for (int j=0; j<=trackstriplist->GetLast(); j++) {
              CandStripHandle * trackstrip = dynamic_cast<CandStripHandle*>(trackstriplist->At(j));

              CandDigitHandleItr digitItr2(trackstrip->GetDaughterIterator());
              CandDigitHandle *digit2 = dynamic_cast<CandDigitHandle*>(digitItr2());

              if(trackstrip->IsEqual(strip)){
                stripIsOnTrack=true;
              }
              if(digit1->GetChannelId()==digit2->GetChannelId()){
                samePixel=true;
              }
            }
          }
        }
        if(!stripIsOnTrack && samePixel){
          dupSlice->RemoveDaughter(strip);
          
        }
        }
      }
      if(dupSlice->GetUidInt()!=slice->GetUidInt()){
        CandTrackSRHandleItr trackItr(cth.GetDaughterIterator());
        for (CandTrackSRHandle *track = trackItr(); track ; track = trackItr()) {
          if(track->GetCandSlice()->GetUidInt()==slice->GetUidInt()){
            CandTrackSRHandle * newtrack = track->DupHandle();
            newtrack->SetCandSlice(dupSlice);
            cth.AddDaughterLink(*newtrack);
            cth.RemoveDaughter(track);
            delete newtrack;
          }
        }
        slicelist->AddDaughterLink(*dupSlice);
        slicelist->RemoveDaughter(slice);
      }
      delete dupSlice;
    } 
  }
}

//----------------------------------------------------------------------------
void   AlgTrackSRList::SpectrometerTracking(TObjArray * track2dlist, 
                                            CandTrackSRListHandle &cth, 
                                            CandSliceHandle *slice,
                                            CandSliceHandle  * dupSlice,
                                            CandContext cx)
{
  MSG("SpectTrackSR",Msg::kDebug) << " In Spectrometer Tracking " << endl;
  // Create Candcontext
  CandContext cxx(this,cx.GetMom());
  
  // Get singleton instance of AlgFactory
  AlgFactory &af = AlgFactory::GetInstance();
  
  CandRecord *candrec = cx.GetCandRecord();
  assert(candrec);
  const VldContext *vldcptr = candrec->GetVldContext();
  assert(vldcptr);
  VldContext vldc = *vldcptr;
  UgliGeomHandle ugh(vldc);

  // generate clones of the digitlist and striplist
  const MomNavigator * mom = cx.GetMom();
  CandRecord *crec = dynamic_cast<CandRecord *>
    (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));

  CandStripListHandle *cslh = dynamic_cast<CandStripListHandle *>
    (crec->FindCandHandle("CandStripListHandle"));
  assert(cslh);
  CandSliceListHandle *cslclh = dynamic_cast<CandSliceListHandle *>
    (crec->FindCandHandle("CandSliceListHandle"));
  assert(cslclh);
  CandDigitListHandle *cdlh = dynamic_cast<CandDigitListHandle *>
    (crec->FindCandHandle("CandDigitListHandle","canddigitlist"));
  assert(cdlh);
  
  if(!fCDLHnew){
    fCDLHnew = cdlh->DupHandle();
    fCDLHnew->SetCandRecord(cdlh->GetCandRecord());
    fCDLHnew->SetName("candspecdigitlist");
  }
  if(!fCSlcLHnew){
    fCSlcLHnew = cslclh->DupHandle();
    fCSlcLHnew->SetCandRecord(cslclh->GetCandRecord());
    fCSlcLHnew->SetName("candspecslicelist");
  }
  if(!fCSLHnew){
    fCSLHnew = cslh->DupHandle();     
    fCSLHnew->SetCandRecord(cslh->GetCandRecord());
    fCSLHnew->SetName("candspecstriplist");
  }

  TIter newdigitItr(fCDLHnew->GetDaughterIterator());
  CandStripHandleItr oldstripItr(cslh->GetDaughterIterator());
  
  Int_t specidir=1;
  //  Int_t minview = min(PlaneView::kU,PlaneView::kV);
  
  oldstripItr.Reset();
  oldstripItr.SetDirection(specidir);
  CandStripHandleItr  newstripItr(fCSLHnew->GetDaughterIterator());     

// generate cluster list for spectrometer hits
// clusters are made for all possible strip end alternatives

//this map stores the correspondence between the new strips we will
// create here and the originals.

  map<CandStripHandle*,CandStripHandle*>newstripmap;
  map<CandStripHandle*,CandStripHandle*>slicemap;

  Int_t specbegplane=1000; 
  TObjArray spectrackclusterlist;

  AlgHandle stripah = af.GetAlgHandle("AlgStripSR","default");
       
  CandStripHandleItr slicestripItr(slice->GetDaughterIterator());
  for (CandStripHandle *strip = slicestripItr(); strip ; strip = slicestripItr()) {
    if (strip->GetPlane()>120 && (strip->GetPlaneView()==PlaneView::kU || strip->GetPlaneView()==PlaneView::kV) ) {
      MSG("SpectTrackSR",Msg::kDebug) << " spect strip " << strip->GetPlane() << " " << strip->GetStrip() << " " << strip->GetCharge() << endl;

      if(strip->GetPlaneView()==PlaneView::kU && strip->GetPlane()>fSliceEndUPlane)
        fSliceEndUPlane=strip->GetPlane();
      if(strip->GetPlaneView()==PlaneView::kV && strip->GetPlane()>fSliceEndVPlane)
        fSliceEndVPlane=strip->GetPlane();

      TIter digitItr(strip->GetDaughterIterator());
      CandDigitHandle *dig=
        dynamic_cast<CandDigitHandle*>(digitItr());           
      const PlexSEIdAltL& altl = dig->GetPlexSEIdAltL();
      int siz = altl.size();
      for (int ind = 0; ind < siz; ++ind) {
        TObjArray diglist;  
        TIter newstripdauItr(strip->GetDaughterIterator());
        while( CandDigitHandle * olddig = dynamic_cast<CandDigitHandle*>(newstripdauItr())){
          CandDigitHandle * newdig = olddig->DupHandle();
          PlexSEIdAltL& newaltl = newdig->GetPlexSEIdAltLWritable(); 
          for (int newind = 0; newind < siz; ++newind) {
            if(newind==ind){
              newaltl[newind].SetWeight((float)1.);
            }
            else{
              newaltl[newind].SetWeight((float)0.);
            }
          }
          newdig->SetCandRecord(olddig->GetCandRecord());
          diglist.Add(newdig);   // diglist does not own newdig.  These must be explicitely deleted 
        }
        cxx.SetDataIn(&diglist);
        CandStripHandle newstrip =
          CandStrip::MakeCandidate(stripah,cxx);
        newstrip.SetCandRecord(cslh->GetCandRecord());
        CandStripHandle *daustrip = new CandStripHandle(newstrip);
        CandStripHandle * nwstrip = dynamic_cast<CandStripHandle *>
          (fCSLHnew->FindDaughter(strip));

        for (Int_t i=0; i<=diglist.GetLast(); i++) {
          CandDigitHandle * deldig =  dynamic_cast<CandDigitHandle*>(diglist.At(i));
          delete deldig;
        }
        
        Bool_t found=false;
        Int_t listlen = spectrackclusterlist.GetLast();
        for (Int_t i=0; i<=listlen; i++) {
          TrackClusterSR* oldcluster = dynamic_cast<TrackClusterSR*>(spectrackclusterlist.At(i));
          assert(oldcluster);
          Int_t minstrip = oldcluster->GetMinStrip();
          Int_t maxstrip = oldcluster->GetMaxStrip();
          if (newstrip.GetPlane()==oldcluster->GetPlane() &&
              (TMath::Abs(minstrip-newstrip.GetStrip())<=fParmTrkClsNSkip+1 || 
               TMath::Abs(maxstrip-newstrip.GetStrip())<=fParmTrkClsNSkip+1)){
            oldcluster->AddStrip(daustrip);
            Int_t istrip=oldcluster->GetNStrip()-1; 
            CandStripHandle * clusterstrip = dynamic_cast<CandStripHandle *>(oldcluster->GetStripList()->At(istrip));
            slicemap[clusterstrip]=strip;
            newstripmap[clusterstrip]=nwstrip;
            found=true;
            break;
          }
        }
        if(!found) {
          MSG("SpectTrackSR",Msg::kDebug) << " spect cluster " << newstrip.GetPlane() << endl;    
          TrackClusterSR *trackcluster = new TrackClusterSR(daustrip, 
                                                            fParmMisalignmentError);
          Int_t istrip=trackcluster->GetNStrip()-1; 
          CandStripHandle * clusterstrip = dynamic_cast<CandStripHandle *>(trackcluster->GetStripList()->At(istrip));
          slicemap[clusterstrip]=strip;
          newstripmap[clusterstrip]=nwstrip;
          spectrackclusterlist.Add(trackcluster);
        
          if (newstrip.GetPlane()<specbegplane)specbegplane = 
                                                       newstrip.GetPlane();
        }
        delete daustrip;
      }
    }
  }   
  fSliceEndPlane = TMath::Max(fSliceEndUPlane,fSliceEndVPlane);

  // now remove all CandStrips in the spectrometer in the slice list
  // and the new strip list.  
       
  slicestripItr.Reset();
  
  MSG("SpectTrackSR",Msg::kDebug) << "# spectrometer track clusters: "
                             << spectrackclusterlist.GetLast() << endl;
  MSG("SpectTrackSR",Msg::kDebug) << "eliminating low pulse height spectrometer clusters" 
                             << endl;
  
 
  TrackClusterSRItr spectclusterItr(&spectrackclusterlist);
  TrackClusterSRKeyFunc *spectclusterKf = spectclusterItr.CreateKeyFunc();
  spectclusterKf->SetFun(TrackClusterSR::KeyFromPlane);
  spectclusterItr.GetSet()->AdoptSortKeyFunc(spectclusterKf);
  spectclusterKf = 0;
  spectclusterItr.SetDirection(specidir);
  spectclusterItr.Reset();
  
  // create vector of hit planes in the spectrometer
  
  vector<Int_t> specplaneindex;
  vector<PlaneView::PlaneView_t> specplaneviewindex;
  Int_t ncount=0;
  Int_t oldipln=0;
  while (TrackClusterSR *tc = spectclusterItr()) {
    Int_t ipln = tc->GetPlane();
    if (!ncount || ipln!=oldipln) {
      specplaneindex.push_back(ipln);
      specplaneviewindex.push_back(tc->GetPlaneView());
    }
    ncount++;
    oldipln = ipln;
  }
  
  // now  add spectrometer track hits to existing tracks
  // no new tracks are formed at this point
  
  Bool_t trackcoverageu(1);
  Bool_t trackcoveragev(1);
  map<Track2DSR*,Int_t> nhitplaneskip;
  
  // loop over hit planes in the spectrometer, from front to back
  
  for (Int_t iplnindx=0; iplnindx<static_cast<Int_t>(specplaneindex.size()) && (trackcoverageu || trackcoveragev);iplnindx++) {
    Int_t ipln = specplaneindex[iplnindx];
    MSG("SpectTrackSR",Msg::kDebug) << "PLANE " << ipln << "/" << specbegplane << "\n";
    if (ipln>=specbegplane) {
      
      Bool_t hasdoubleended(0);
      
      // loop over spectrometer clusters, finding those in plane ipln
      // these clusters are held in tclusterpln
      
      TObjArray tclusterpln;
      for (Int_t i=0; i<=spectrackclusterlist.GetLast(); i++) {
        TrackClusterSR* tc = dynamic_cast<TrackClusterSR*>(spectrackclusterlist.At(i));
        if (!hasdoubleended && tc->IsDoubleEnded()) {
          hasdoubleended = 1;
        }
        
        if (tc->GetPlane()==ipln) {
          tclusterpln.Add(tc);
          if (tc->GetPlaneView()==PlaneView::kU) {
            trackcoverageu = 0;
          }
          if (tc->GetPlaneView()==PlaneView::kV) {
            trackcoveragev = 0;
          }
        }
      }
      Int_t useph = 0;
      if (fParmIsCosmic==1 && hasdoubleended) {
        useph = 1;
      }
      // loop over tracks. Each track will be extrapolated to and compared
      // with the position of each cluster on plane ipln
      
      for (Int_t itrk=0; itrk<=track2dlist->GetLast(); itrk++) {
        Track2DSR *thistrack = dynamic_cast<Track2DSR*>(track2dlist->At(itrk));
        MSG("SpectTrackSR",Msg::kDebug) << "track " << itrk << "/" << track2dlist->GetLast() << "\n";
        Double_t residParams[] = {1.e6,1.e6};
        TrackClusterSR *besttc=0;
        PlaneView::PlaneView_t tcplaneview = PlaneView::kUnknown;
        Bool_t shouldhit(kFALSE); // true is same view as track
        
        // loop over track end interval to extrapolate from
        //  Int_t ibefpln=0;
        for (Int_t ibefpln=0; ibefpln<=fParmTrk2DNSkipRemove && ibefpln<thistrack->GetLast() && !besttc; ibefpln++) {
          Int_t numSkippedHit = 0;
          TrackClusterSR *tct = 
            (thistrack->GetCluster(ibefpln));
          
          Int_t slopelookup=ibefpln;
          if(slopelookup<thistrack->GetLast())slopelookup++;
          Double_t trackSlope = thistrack->GetForwardSlope(ibefpln+1);
          Int_t deltaPlane=0;
          Int_t dview=0;
          Int_t numSkippedActive=0; // number of active planes of same view
          
          // loop over all clusters in plane ipln 
          
          for (Int_t itc=0; itc<=tclusterpln.GetLast(); itc++) {
            TrackClusterSR *tc = dynamic_cast<TrackClusterSR*>(tclusterpln.At(itc));
            
            if(!itc){
              deltaPlane = tc->GetPlane()-tct->GetPlane();
              tcplaneview = tc->GetPlaneView();
              dview = tcplaneview-tct->GetPlaneView();
              
              if (dview==0) shouldhit = kTRUE;
              if (shouldhit && deltaPlane>0) {
                Float_t upos=0;
                Float_t vpos=0;
                Float_t uslope=0;
                Float_t vslope=0;
                if(tcplaneview==PlaneView::kU){
                  upos = tct->GetTPos();
                  vpos = fHoughVIntercept + tct->GetZPos()*fHoughVSlope;
                  uslope=trackSlope;
                  vslope=fHoughVSlope;
                }
                else{
                  vpos = tct->GetTPos();
                  upos = fHoughUIntercept + tct->GetZPos()*fHoughUSlope;
                  vslope=trackSlope;
                  uslope=fHoughUSlope;
                }
                numSkippedActive=FindNumSkippedPlanes(tc->GetPlane(),
                                                      tct->GetPlane(),
                                                      upos,vpos,
                                                      tct->GetZPos(),
                                                      uslope,
                                                      vslope,
                                                      specidir,
                                                      tct->GetPlaneView());
                
              }
            }
            Int_t nhit = thistrack->GetLast()-ibefpln+1;
            Float_t aveSpacing=2;
            if(thistrack->GetLast()>0){
              aveSpacing = TMath::Abs(thistrack->GetEndPlane()-thistrack->GetBegPlane())/ (thistrack->GetLast());
            }         
            MSG("SpectTrackSR",Msg::kDebug) << "nhit, nplaneskip = " << nhit << "/" << fParmTrk2DNContiguous  
                                            << "  " << numSkippedActive << "/" << fParmTrk2DNSkip 
                                            << " dview/dplane " << dview << "/" << deltaPlane 
                                            << " track plane/min track plane " << tct->GetPlane() 
                                            << "/" << 120 - aveSpacing*fParmTrk2DNSkip*2 << endl;
            
            if (dview==0 && 
                deltaPlane>=0 &&  
                numSkippedActive<=fParmTrk2DNSkip && 
                //  tct->GetPlane()>=(120-aveSpacing*fParmTrk2DNSkip*2) && 
                (nhit>=fParmTrk2DNContiguous || numSkippedActive==0)) { 
              
              MSG("SpectTrackSR",Msg::kDebug) << " checking SE alternative: Tpos: " << tc->GetTPos() <<  endl;
              if (tcplaneview==PlaneView::kU) {
                trackcoverageu = 1;
                MSG("SpectTrackSR",Msg::kDebug) << "clearing U trackcoverage\n";
              }
              if (tcplaneview==PlaneView::kV) {
                trackcoveragev = 1;
                MSG("SpectTrackSR",Msg::kDebug) << "clearing V trackcoverage\n";
              }
              
              //reset the residParams because you are on a new plane
              Bool_t isSpectrometerPlane=true;
              
              if(IsBestClusterInPlane(tc,tct, numSkippedHit,
                                      nhit, ibefpln, trackSlope,thistrack, 
                                      
                                      residParams, 
                                      deltaPlane,isSpectrometerPlane))  besttc = tc;
            }
          }
                } // end ibefpln
        
        if (besttc) { 
          TrackClusterSR * newtrackcluster = 0;
          
          // loop over strips in this cluster, and for each....
          //    loop over digit daughters, and find on new digit list
          //     3 possibilities exist...
          //     new digit has zero weights - clone with correct weights
          //         and delete zero weight version.
          //     new digit has weights already set (correctly) do nothing
          //     new digit has weight set, but pointing to a different strip
          //         in this case, make new CandDigit.  
          
          for (Int_t istrip=0; istrip<besttc->GetNStrip(); istrip++) {
            CandStripHandle * strip = dynamic_cast<CandStripHandle *>(besttc->GetStripList()->At(istrip));
            assert(strip);
            
            // now build a new strip to be added to strip list
            
            TObjArray cdhAr;
            cdhAr.Clear();
            Bool_t deleteold=true;
            // loop over digit daughters of strip
            
            TIter digitItr(strip->GetDaughterIterator());
            while (CandDigitHandle *dig=dynamic_cast<CandDigitHandle*>(digitItr())) {    
              newdigitItr.Reset();
              Bool_t moddig=false;
              Bool_t prevadded=false;
              CandDigitHandle * founddig=0;
              
              // find this digit in the new digit list
              
              while (CandDigitHandle *newdig=dynamic_cast<CandDigitHandle*>(newdigitItr())){
                if(dig->GetChannelId()==newdig->GetChannelId() && dig->GetRawDigitIndex()==newdig->GetRawDigitIndex()){ 
                  founddig=newdig;                     
                  const PlexSEIdAltL& oldaltl = newdig->GetPlexSEIdAltL();
                  
                  // digit is found - check weights
                  
                  if(oldaltl.GetBestWeight()==0) {
                    // if best weight is zero, we make a version of this
                    // digit with correct weights first duplicating the handle.
                    CandDigitHandle * dupdig = founddig->DupHandle(); 
                    PlexSEIdAltL& dupaltl = dupdig->GetPlexSEIdAltLWritable();
                    moddig=true;
                    int siz = dupaltl.size();
                    for (int ind = 0; ind < siz; ++ind) {
                      if(dupaltl[ind].GetSEId().GetStrip()==
                         strip->GetStrip()){
                        dupaltl[ind].SetWeight((float)1.);
                      }
                      else{
                        dupaltl[ind].SetWeight((float)0.);
                      }
                    }
                    
                    dupdig->SetCandRecord(founddig->GetCandRecord());
                    fCDLHnew->AddDaughterLink(*dupdig);
                    CandDigitHandle * daudig = dynamic_cast<CandDigitHandle *>(fCDLHnew->FindDaughter(dupdig));
                    assert(daudig);
                    fCDLHnew->RemoveDaughter(newdig);
                    cdhAr.Add(daudig);
                    delete dupdig;
                  }
                  else {
                    
                    // this digit has already been set to the correct SEId, 
                    // so just leave it alone.
                    
                    if(oldaltl.GetBestWeight()==1.0 && oldaltl.GetBestItem().GetSEId().GetStrip()==strip->GetStrip()){
                      prevadded=true;
                      cdhAr.Add(newdig);
                    }
                  }
                  break;
                }
              }
              
              if(founddig && !moddig && !prevadded){ 
                // need to construct new CandDigit, as there is
                // no existing digit which has proper altl weights.  Since
                // the weights have apparently already been set in the
                // found digit (probably associated 
                // with an earlier track),
                // leave it in place.
                deleteold=false;
                CandDigitHandle * dupdig = founddig->DupHandle(); 
                PlexSEIdAltL& dupaltl = dupdig->GetPlexSEIdAltLWritable();
                int siz = dupaltl.size();
                for (int ind = 0; ind < siz; ++ind) {
                  if(dupaltl[ind].GetSEId().GetStrip()==
                     strip->GetStrip()){
                    dupaltl[ind].SetWeight((float)1.);
                  }
                  else{
                    dupaltl[ind].SetWeight((float)0.);
                  }
                }
                dupdig->SetCandRecord(founddig->GetCandRecord());
                fCDLHnew->AddDaughterLink(*dupdig);
                CandDigitHandle * daudig = dynamic_cast<CandDigitHandle *>(fCDLHnew->FindDaughter(dupdig));
                if(daudig)cdhAr.Add(daudig);
                delete dupdig;
              }       
            }
            
            // we now construct the strip, and add to striplist and slice
            // daughter list
            
            AlgHandle stripah = af.GetAlgHandle("AlgStripSR","default");
            if(cdhAr.GetLast()>=0){
              cxx.SetDataIn(&cdhAr);
              CandStripHandle striphandle =
                CandStrip::MakeCandidate(stripah,cxx);
              striphandle.SetCandRecord(cslh->GetCandRecord());

              fCSLHnew->AddDaughterLink(striphandle);           
              dupSlice->AddDaughterLink(striphandle);
              CandStripHandle * daustrip = dynamic_cast<CandStripHandle *>(fCSLHnew->FindDaughter(&striphandle));

              assert (daustrip);
              // add to track cluster list
              if(istrip==0 || !newtrackcluster){
                newtrackcluster = new TrackClusterSR(daustrip, 
                                                     fParmMisalignmentError);
              }
              else{
                newtrackcluster->AddStrip(daustrip);
              }

            }
          }
          if(newtrackcluster){
            // add this cluster to track
            TrackClusterSR *tc0 = thistrack->GetCluster(0);
            Double_t slope = (newtrackcluster->GetTPos()-tc0->GetTPos())/
              (newtrackcluster->GetZPos()-tc0->GetZPos());
            // add the full spectrometer cluster to the track, and to the
            // track cluster list
            thistrack->Add(newtrackcluster,slope);
            cth.AddTrackCluster(newtrackcluster);
            delete newtrackcluster;
          }
          nhitplaneskip[thistrack]=0;
        }
        else if (shouldhit && ipln>thistrack->GetCluster(thistrack->GetLast())->GetPlane()) {
          nhitplaneskip[thistrack]++;
        }
      }
    }
  }
  spectrackclusterlist.Delete();

  if(dupSlice->GetUidInt()!=slice->GetUidInt()){
    fCSlcLHnew->AddDaughterLink(*dupSlice);
    fCSlcLHnew->RemoveDaughter(slice);
  }
  // check each track for outlier hit at two track ends, and remove if found
  for (Int_t itrk=0; itrk<=track2dlist->GetLast(); itrk++) {
    Track2DSR *track = dynamic_cast<Track2DSR*>(track2dlist->At(itrk)); 
    if(track->GetLast()>2){
      Float_t gap1=TMath::Abs(track->GetCluster(0)->GetPlane()-track->GetCluster(1)->GetPlane());
      Float_t prevgap1=TMath::Abs(track->GetCluster(2)->GetPlane()-track->GetCluster(1)->GetPlane());
      gap1=gap1/prevgap1;
      Float_t gap2=TMath::Abs(track->GetCluster(track->GetLast())->GetPlane()-track->GetCluster(track->GetLast()-1)->GetPlane());
      Float_t prevgap2=TMath::Abs(track->GetCluster(track->GetLast()-1)->GetPlane()-track->GetCluster(track->GetLast()-2)->GetPlane());
      gap2=gap2/prevgap2;
      Int_t gapfactor=2;
      if(gap1>fParmTrk2DNSkipHit*gapfactor)track->RemoveAt(0);
      if(gap2>fParmTrk2DNSkipHit*gapfactor)track->RemoveAt(track->GetLast());
      track->Compress();
    }   
  }
}
//----------------------------------------------------------------------
//this method determines the direction the track is going based on 
//timing.  it returns a 1 if the direction goes in +z and -1 if in -z
Int_t AlgTrackSRList::FindTimingDirection(CandStripHandleItr &allStripItr, 
                                          Float_t uTrackSlope, 
                                          Float_t uTrackIntercept,
                                          Float_t vTrackSlope, 
                                          Float_t vTrackIntercept)
{
  UgliGeomHandle ugh(fVldc);

  Int_t idir = 1;

  ArrivalTime arrtime;

  // determine directionality using timing - this should really
  //only be needed for cosmic ray events, as beam events always
  //increase in z with increasing time.
  
  //we could try to determine the timing for each side (E, W) of each
  //electronics crate.  that would be the most accurate way, but it 
  //would also take a lot of time and iterations.  instead, we can 
  //probably get away with just doing a global fit to all datapoints
  //since we do have timing constants that take care of the side to side
  //and crate to crate offsets.  even if a bit of timing hardware is 
  //exchanged, this wont really affect our ability to go back and figure
  //out the right calibrations later.
  
  //make some arrays to do the fit. allow a maximum of 1000 points
  //to be included in the fit.  there may be events with more than
  //1000 digits, but that should be enough to tell you which way the
  //event is going.
  Double_t zpos[1000];
  Double_t time[1000];
  Double_t weight[1000];
  
  for (int i=0; i<1000; i++) {
    zpos[i] = 0.;
    time[i] = 0.;
    weight[i] = 0.;
  }
  
  //declare variables used in the while loop below
  //distFromCenter is the offset in u or v from the center of the strip
  //halflength is the strip halflength
  //upos and vpos are obvious
  //fiberDist is the total distance traveled in fiber 
  Int_t digitCtr = 0;
  Double_t distFromCenter = 0.;
  Double_t halflength = 0.;
  Double_t upos = 0.;
  Double_t vpos = 0.;
  Double_t fiberDist = 0.;
  CandStripHandle *strip = 0;
  CandDigitHandle *digit = 0;

  PlaneView::PlaneView_t view0 = PlaneView::kU;
  PlaneView::PlaneView_t view1 = PlaneView::kV;

  //minimum weight for a signal in the timing fit
  // maximum weight corresponds to 10 p.e.; reason for this is to minimize
  // weight of showers in which transverse spread can have negative effect
  // on timing
  const Double_t min_wgt = 0.4;
  
  allStripItr.ResetFirst();
  MSG("AlgTrackSRList", Msg::kDebug) << "-2 fit over " << digitCtr << " points" << endl;

  while( (strip = allStripItr()) && digitCtr<1000){
    
    //get the Plex and Ugli information for this strip to be used when 
    //figuring out the fiber lengths
    PlexStripEndId stripid = strip->GetStripEndId();
    UgliScintPlnHandle planehandle = ugh.GetScintPlnHandle(stripid);
    TIter digitItr(strip->GetDaughterIterator());
    UgliStripHandle striphandle = ugh.GetStripHandle(stripid);
    halflength = striphandle.GetHalfLength();
    upos = uTrackIntercept+uTrackSlope*strip->GetZPos();
    vpos = vTrackIntercept+vTrackSlope*strip->GetZPos();
    
    //loop over the digits associated with this strip and only use those
    //which dont have the demux veto flag set - those that do were not
    //used in the demuxing solution
    for (int j=0; j<strip->GetNDigit() && digitCtr<1000; ++j) {
      digit = dynamic_cast<CandDigitHandle*>(digitItr());
      
//------------>>>>coud also require a minimum pulseheight to be included in the fit
      if (!digit->GetPlexSEIdAltL().GetDemuxVetoFlag()) {
        //get the z position of the plane
        zpos[digitCtr] = (Double_t)planehandle.GetZ0();
        
        PlexSEIdAltL altl = digit->GetPlexSEIdAltL();
        StripEnd::StripEnd_t stripEnd = altl.GetEnd();
        distFromCenter = 0.;
        
        // if we're in a U plane we need to account for the V offset
        // and vice versa
        if (digit->GetPlexSEIdAltL().GetPlaneView() == view1) {
          distFromCenter = 
            (digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kNegative) ?  upos : -upos;
        }
        else if (digit->GetPlexSEIdAltL().GetPlaneView() == view0) {
          distFromCenter = 
            (digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kNegative) ? -vpos :  vpos;
        }
        
        fiberDist = ((halflength + distFromCenter) +
                     striphandle.ClearFiber(stripEnd) + 
                     striphandle.WlsPigtail(stripEnd));
        
        time[digitCtr] = (strip->GetTime(digit->GetPlexSEIdAltL().GetEnd()) 
                          - (fiberDist/PropagationVelocity::Velocity(fSimFlag)));
        
        weight[digitCtr] = min(arrtime.Weight(digit->GetCharge(CalDigitType::kPE)),min_wgt);
        ++digitCtr;
        
        MSG("AlgTrackSRList", Msg::kDebug) << "-1 fit over " << digitCtr << " points" << endl;

      }//end if this is a non-vetoed digit
    }//end loop over digits
  }//end loop over strips to fill timing arrays
  
  //reset the iterator over all strips
  allStripItr.ResetFirst();
  
  //arrays to hold the slope, intercept, and chi^2 of the 
  //least squares fit for the timing
  Double_t parm[2],eparm[2];
  
  //set the maximum residual to the input parameter+1ns to 
  //get ready for the while loop below
  Double_t maxresid = fParmMaxTimingResid+1.*Munits::ns;
  Double_t resid = 0.;
  Int_t imaxresid = 0;
  Int_t nremove=0;
  Int_t fitflag=0;
  
  MSG("AlgTrackSRList", Msg::kDebug) << "0 fit over " << digitCtr << " points" << endl;
  
  while(digitCtr-nremove>4 && maxresid>fParmMaxTimingResid && !fitflag){
    MSG("AlgTrackSRList", Msg::kDebug) << "1 fit over " << digitCtr << " points" << endl;
    fitflag = LinearFit::Weighted(digitCtr,zpos,time,weight,parm,eparm);
    maxresid = 0.;
    imaxresid = 0;
    for (int i=0; i<digitCtr; ++i) {
      resid = parm[0]+parm[1]*zpos[i]-time[i];
      if (weight[i]>0. && TMath::Abs(resid)>maxresid) {
        maxresid = TMath::Abs(resid);
        imaxresid = i;
      }
    }
    if (maxresid>fParmMaxTimingResid) {
      weight[imaxresid] = 0.;
      ++nremove;
    }
    MSG("AlgTrackSRList", Msg::kDebug) << "2 fit over " << digitCtr << " points" << endl;
  }//end loop to fit the timing values

  MSG("AlgTrackSRList", Msg::kVerbose) << "timing slope is " << parm[1] << endl;

  if(parm[1]<0.){
    idir=-1;
  }

  return idir;
}

//----------------------------------------------------------------------
//method to find the direction of the event wrt z
Int_t AlgTrackSRList::FindTimingDirection(TrackClusterSRItr &clusterItr)
{
  Int_t idir = 1;

  ArrivalTime arrtime;

  // determine directionality using timing - this should really
  //only be needed for cosmic ray events, as beam events always
  //increase in z with increasing time.
  
  //we could try to determine the timing for each side (E, W) of each
  //electronics crate.  that would be the most accurate way, but it 
  //would also take a lot of time and iterations.  instead, we can 
  //probably get away with just doing a global fit to all datapoints
  //since we do have timing constants that take care of the side to side
  //and crate to crate offsets.  even if a bit of timing hardware is 
  //exchanged, this wont really affect our ability to go back and figure
  //out the right calibrations later.
  
  //make some arrays to do the fit. allow a maximum of 1000 points
  //to be included in the fit.  there may be events with more than
  //1000 digits, but that should be enough to tell you which way the
  //event is going.
  Double_t zpos[1000];
  Double_t time[1000];
  Double_t weight[1000];
  
  for (int i=0; i<1000; i++) {
    zpos[i] = 0.;
    time[i] = 0.;
    weight[i] = 0.;
  }
  
  //minimum weight for a signal in the timing fit
  // maximum weight corresponds to 10 p.e.; reason for this is to minimize
  // weight of showers in which transverse spread can have negative effect
  // on timing
  const Double_t min_wgt = 0.4;
  TrackClusterSR *cluster = 0;

  Int_t clusterCtr = 0;

  clusterItr.ResetFirst();
  while( (cluster = clusterItr()) && clusterCtr<1000){
    
    //get the z position of the plane
    zpos[clusterCtr] = (Double_t)cluster->GetZPos();
        
    time[clusterCtr] = cluster->GetBegTime();
                
    weight[clusterCtr] = min(arrtime.Weight(cluster->GetCharge()),min_wgt);
    ++clusterCtr;

  }//end loop over clusters to fill timing arrays
  
  //reset the iterator over all strips
  clusterItr.ResetFirst();
  
  //arrays to hold the slope, intercept, and chi^2 of the 
  //least squares fit for the timing
  Double_t parm[2],eparm[2];
  parm[0]=0;
  parm[1]=0;
  eparm[0]=0;
  eparm[1]=0;
  //set the maximum residual to the input parameter+1ns to 
  //get ready for the while loop below
  Double_t maxresid = fParmMaxTimingResid+1.*Munits::ns;
  Double_t resid = 0.;
  Int_t imaxresid = 0;
  Int_t nremove=0;
  Int_t fitflag = 0;

  while(clusterCtr-nremove>4 && maxresid>fParmMaxTimingResid && !fitflag){
    fitflag = LinearFit::Weighted(clusterCtr,zpos,time,weight,parm,eparm);
    maxresid = 0.;
    imaxresid = 0;
    for (int i=0; i<clusterCtr; ++i) {
      resid = parm[0]+parm[1]*zpos[i]-time[i];
      if (weight[i]>0. && TMath::Abs(resid)>maxresid) {
        maxresid = TMath::Abs(resid);
        imaxresid = i;
      }
    }
    if (maxresid>fParmMaxTimingResid) {
      weight[imaxresid] = 0.;
      ++nremove;
    }
  }//end loop to fit the timing values

  MSG("AlgTrackSRList", Msg::kDebug) << "cluster timing slope is " << parm[1] << endl;

  if(parm[1]<0.){
    idir=-1;
  }

  return idir;
}

//-----------------------------------------------------------------------
//method to identify seed clusters and put them into a trackList as the
//first clusters in a track, returns the number of tracks found

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

  TrackClusterSR *cluster = 0;
  TrackClusterSR *nextCluster=0;
  TrackClusterSRItr nextclusterItr2(clusterItr);

  Int_t prevTrackCount = trackList->GetLast()+1;
  MSG("TrackSR", Msg::kDebug) << "previously " << prevTrackCount 
                             << " tracks in list" << endl;
 //set the direction of the iterator based on the timing direction.
  //then the first plane represent by clusters is your first plane
  //that is a possible seed hit.
  if(direction == 1) clusterItr.ResetFirst();
  else if(direction == -1) clusterItr.ResetLast();

  vector<Int_t> planeindex;
  Int_t ncount=0;
  Int_t oldipln=0;
  while (TrackClusterSR *tc = clusterItr()) {
    Int_t ipln = tc->GetPlane();

    if(tc->IsValid()){
      if (!ncount || ipln!=oldipln ) {
        planeindex.push_back(ipln);
        MSG("TrackSR", Msg::kDebug) << " valid cluster on plane " << ipln << endl;
      }
      ncount++;
      oldipln = ipln;
    }
  }
  
  //difference in plane number between plane of current cluster
  //and one being looked at
  Int_t deltaPlane = 0;

  //difference in time between current cluster and one being looked at
  Double_t deltaTime = 0.;
  
  //difference in z pos and transverse position between current cluster
  //and one being looked at
  Double_t deltaTPos = 0.;
  Double_t deltaZPos = 0.;

  //keep track of which valid cluster you are on
  Int_t validCtr = 0;
  //loop over all the clusters in the list.  only those that are still
  //valid will be considered to be used as new seed clusters.  non-valid
  //clusters have alread been used in a track
  for (Int_t iplnindx=0; iplnindx<static_cast<Int_t>(planeindex.size());iplnindx++) {
    Int_t ipln = planeindex[iplnindx];
    MSG("TrackSR",Msg::kDebug) << "PLANE " << ipln  << "\n";
    //reset the direction of the iterator.
    if(direction == 1) clusterItr.ResetFirst();
    else if(direction == -1) clusterItr.ResetLast();
    while( (cluster = clusterItr()) ){
      if(cluster->GetPlane()==ipln && cluster->IsValid()){
        //if you have a valid cluster, it could be a seed hit
         
        MSG("TrackSR",Msg::kDebug) << cluster->GetPlane() 
                                   << " " << cluster->GetTPos() << " " 
                                   << cluster->GetZPos()
                                   << " " << cluster->GetCharge() << " " 
                                   << (Int_t)cluster->IsValid() << " " 
                                   << cluster->GetNStrip() << endl;
        
        //seed clusters are those which are well separated from clusters
        //preceding them in the list in terms of plane number.  they must 
        //also be well separated in time and transverse position
        
        if(iplnindx==0){
          ++validCtr;
          //first valid plane in the bunch, call it begPlane
          begPlane = cluster->GetPlane();
          //first valid cluster in the view, it is a seed cluster
          Track2DSR *track = new Track2DSR();
          MSG("TrackSR",Msg::kDebug) << " Found Seed hit " << cluster->GetNStrip() << "/" << fParmSingleHitDef ;
          //set the fIterate variable for track so it can id
          //the iteration it was created in
          track->fIterate = iterate;
          track->SetDirection(direction);
          if(!cluster->IsWide()){
            MSG("TrackSR",Msg::kDebug) << " OK as seed" << endl;
            seedClusters->AddLast(cluster);
          }
          //too many hits in the cluster to use as a seed cluster
          //set the cluster as non-valid so it doesnt get used
          //in the initial finding somehow.  later when we loop over
          //all tracks to fill in the gaps it could be used
          else{
            MSG("TrackSR",Msg::kDebug) << " too wide for seed " << endl;
            track->fIterate = -1;
            cluster->SetIsValid(false);
          }
          track->Add(cluster,houghSlope);
          track->SetHoughSlope(houghSlope);
          track->SetHoughIntercept(houghIntercept);
          trackList->AddLast(track);
        }
        else{
          Bool_t isolated=true;
    //reset the direction of the iterator.
          if(direction == 1) nextclusterItr2.ResetFirst();
          else if(direction == -1) nextclusterItr2.ResetLast();
          while( (nextCluster = nextclusterItr2()) ){
            if(nextCluster->GetPlane()==begPlane && nextCluster->IsValid()){
              deltaPlane = TMath::Abs(cluster->GetPlane() - nextCluster->GetPlane());
              deltaTime = TMath::Abs(cluster->GetBegTime() - nextCluster->GetBegTime());
              deltaZPos = TMath::Abs(cluster->GetZPos() - nextCluster->GetZPos());
              deltaTPos = min(TMath::Abs(cluster->GetTPos() - nextCluster->GetTPos() 
                                     + deltaZPos*houghSlope),TMath::Abs(cluster->GetTPos() - nextCluster->GetTPos()));
              
              MSG("TrackSR", Msg::kDebug) << "current cluster plane = " 
                
                                          << cluster->GetPlane()
                                          << " prev cluster plane = " 
                                          << nextCluster->GetPlane()
                                          << endl << " deltaPlane = " << deltaPlane 
                                          << " / " << fParmTrk2DPlnEnd << endl  
                                          << " deltaTime = " << deltaTime << "/ " 
                                          << fParmHitNTime << "/" << fParmHitTime
                                          << " deltaTPos = " << deltaTPos << " / " 
                                          << fParmTrk2DWin0*deltaZPos << endl;
              
              if(deltaPlane<=fParmTrk2DPlnEnd 
                 && deltaTime>fParmHitNTime && deltaTime<fParmHitTime
                 && deltaTPos<fParmTrk2DWin0*deltaZPos)
                {
                  isolated=false;
                  break;
                }
            }
          }
          if(isolated){
            //make a new track
            Track2DSR *track = new Track2DSR();
            MSG("TrackSR",Msg::kDebug) << " Found Seed hit " << endl;
            //set the fIterate variable for track so it can id
            //the iteration it was created in
            track->fIterate = iterate;
            track->SetDirection(direction);
            if(!cluster->IsWide()){
              seedClusters->AddLast(cluster);
            }
            //too many hits in the cluster to use as a seed cluster
            //set the cluster as non-valid so it doesnt get used
            //in the initial finding somehow.  later when we loop over
            //all tracks to fill in the gaps it could be used
            else{
              track->fIterate = -1;
              cluster->SetIsValid(false);
            }
            track->Add(cluster,houghSlope);
            track->SetHoughSlope(houghSlope);
            track->SetHoughIntercept(houghIntercept);
            trackList->AddLast(track);
          }       
        }
      }
    }
    begPlane=ipln;
  }

  //loop over the tracks one more time to get rid of the flagged tracks
  //i could have maybe did this in the above loop, but its not clear
  //that i wouldnt end up deleting an object before the last time it would
  //be accessed by the loop

  Track2DSR *track1 = 0;
  for(Int_t i = 0; i<=trackList->GetLast(); ++i){
    track1 = dynamic_cast<Track2DSR *>(trackList->At(i));
    if(track1->fIterate < 0){
      trackList->RemoveAt(i);
      delete track1;
    }
  }//end loop to remove tracks
  
  trackList->Compress();
  return trackList->GetLast()-prevTrackCount;
}

//------------------------------------------------------------------------
//this method adds clusters to the input tracks.  before calling this method
//we have ensured that it was created in the current iteration of looking
//for and forming tracks, and we of course are only adding clusters to 
//the track which come from the same view as the track
Int_t AlgTrackSRList::AddClustersToTracks(Track2DSRItr &trackItr, 
                                          TrackClusterSRItr &planeClusterItr,
                                          Int_t direction, Int_t iterate,
                                          std::vector<Int_t>& hitPlanes, Int_t trk2dmin)
{

  trackItr.ResetFirst();

  Int_t currentPlane = 0;
  Int_t hitPlanesSkippedInTrack = 0;
  TrackClusterSR *besttc=0;  //place holder of the cluster to add
  TrackClusterSR *prevCluster=0;  //previous cluster in track to current plane
  TrackClusterSR *planeCluster=0; //cluster of current plane we are looking at
  Track2DSR *track = 0; //current track

  //set the direction of your iterators
  if(direction == 1) planeClusterItr.ResetFirst();
  else if(direction == -1) planeClusterItr.ResetLast();
 
  Int_t planesBefore = 0;

  //numSkippedActive = # of active planes we cant ignore between
  //the current plane and the plane of the last cluster in the current 
  //track.  numSkippedHit = # of hit planes in this event skipped
  Int_t numSkippedActive = 0;
  Int_t numSkippedHit = 0;

  //this is the slope of the track at the cluster you are looking at
  Double_t trackSlope = 0.;

  //hold the values of the fit for the clusters in the following array
  //residParams[0] is the best residual without using signal weights, 
  //residParams[1] is the best residual using signal weights
  Double_t residParams[] = {1.e6,1.e6};

  //number of planes between current cluster and previous one already in track
  Int_t deltaPlane = 0;
  Int_t trackCtr = 0;
  
  //# of planes between the one with the last cluster added to the track
  //and the plane in the track holding prevCluster
  Int_t nHit = 0;

  //make an array to hold the index of hitPlanes for a given plane
  Int_t planeIndex[500];
  for(Int_t i = 0; i < 500; ++i){
    planeIndex[i] = 0;
  }
  for(Int_t ii = 0; ii < static_cast<Int_t>(hitPlanes.size()); ++ii){
   planeIndex[hitPlanes[ii]] = ii;
  
  }
  Int_t index = 0;

  //loop over the tracks in the list - this way you add clusters from
  //the current plane to each track and save cpu time vs looping over
  //each track and then every plane in the list
  trackCtr = 0;
  while( (track = trackItr()) ){   
    MSG("TrackSR", Msg::kDebug) << "iterate = " << iterate 
                                << " track->fIterate = " << track->fIterate 
                                << " track # " << trackCtr 
                                << " track beg plane = " << track->GetBegPlane() 
                                << " track is bad = " << (Int_t)track->IsBad() 
                                << endl;
    index = planeIndex[track->GetBegPlane()];  
     
    //if this track was created in this iteration and was not
    //previously marked as a bad track
    if( track->fIterate == iterate && !track->IsBad() ){
      
      MSG("AlgTrackSRList", Msg::kVerbose) << "on a good track in iteration " 
                                         << iterate << endl;
      planesBefore = 0;
      planeClusterItr.ResetFirst();

      //loop over the planes hit in this slice
      for(Int_t plnIndex=index; plnIndex<static_cast<Int_t>(hitPlanes.size()); 
          ++plnIndex){
        
        currentPlane = hitPlanes[plnIndex];
        deltaPlane = 0;
        MSG("TrackSR", Msg::kDebug) << "currentPlane = " << currentPlane << endl;                       
        numSkippedActive = 0;
        numSkippedHit = 0;
        
        //have to reset planesBefore for each new track
        planesBefore = 0;
                
        //reset besttc because you are on a new plane
        besttc = 0;
        
        // make sure that the current plane is downstream of 
        // the first cluster in the list of clusters for this time through the finder

        if(currentPlane*direction
           >(direction*track->GetCluster(track->GetLast())->GetPlane())){ 
              
          //loop from the current last cluster in the track to the beginning of the
          //track if necessary to see if the cluster of the next plane matches
          //with any of them

          MSG("TrackSR", Msg::kDebug) << "planesBefore = " 
                                      << planesBefore << endl;

          //loop backwards over the clusters in the track
          while( planesBefore<= fParmTrk2DNSkipRemove 
                 && track->GetLast()-planesBefore>=0
                 && !besttc){

            //get the cluster you are interested in
            prevCluster = track->GetCluster(track->GetLast()-planesBefore);
            
            MSG("TrackSR", Msg::kDebug) << track->GetForwardSlope(track->GetLast()-planesBefore)
                                        << " " << planesBefore << " " 
                                        << track->GetLast()+1 << endl;

            trackSlope = track->GetForwardSlope(track->GetLast()-planesBefore);
            Float_t upos=0;
            Float_t vpos=0;
            Float_t uslope=0;
            Float_t vslope=0;
            if(prevCluster->GetPlaneView()==PlaneView::kU){
              upos = prevCluster->GetTPos();
              vpos = fHoughVIntercept + prevCluster->GetZPos()*fHoughVSlope;
              uslope=track->GetForwardSlope(track->GetLast()-planesBefore);
              vslope=fHoughVSlope;
            }
            else{
              vpos = prevCluster->GetTPos();
              upos = fHoughUIntercept + prevCluster->GetZPos()*fHoughUSlope;
              vslope=track->GetForwardSlope(track->GetLast()-planesBefore);
              uslope=fHoughUSlope;
            }
            numSkippedActive=FindNumSkippedPlanes(currentPlane,
                                                  prevCluster->GetPlane(),
                                                  upos,vpos,
                                                  prevCluster->GetZPos(),
                                                  uslope,
                                                  vslope,
                                                  direction,
                                                  prevCluster->GetPlaneView());
            
            deltaPlane = TMath::Abs(prevCluster->GetPlane()-currentPlane);
            
            MSG("TrackSR", Msg::kDebug) <<"current plane " 
                                        << currentPlane 
                                        <<" prev plane " 
                                        << prevCluster->GetPlane() 
                                        <<" active skipped = " 
                                        << numSkippedActive
                                        <<" deltaPlane = " 
                                        << deltaPlane << endl;
            
            //the +1 below accounts for the fact that get last
            //returns a number starting the count in the array at 0
            nHit = track->GetLast()-planesBefore+1;
            
            //if the distance between the plane of the previous cluster and
            //the current plane seems reasonable, then look for clusters 
            //to add from the current plane.  make sure you arent looking 
            //at clusters from the same plane though
            if(numSkippedActive<=fParmTrk2DNSkip
               && (nHit>=fParmTrk2DNContiguous || numSkippedActive==0)
               && currentPlane != prevCluster->GetPlane() ){
              
              //reset the residParams because you are on a new plane
              residParams[0] = 1.e6;
              residParams[1] = 1.e6;

             
              //slice the planeClusterItr to get only those clusters on the current plane
              planeClusterItr.GetSet()->Slice();
              planeClusterItr.GetSet()->Slice(currentPlane);     
              //loop over the clusters in the current plane
              while( planeClusterItr.IsValid() ){
                planeCluster = planeClusterItr.Ptr();
                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;
                                
                Bool_t isSpectrometerPlane = false; 
                if(fDetector==Detector::kNear 
                   && (planeCluster->GetPlane()>=121)){
                  isSpectrometerPlane=true;
                }
                if(planeCluster->IsValid()){
                  if(IsBestClusterInPlane(planeCluster,prevCluster, 
                                          numSkippedHit,
                                          nHit, planesBefore,trackSlope, track, residParams, 
                                          deltaPlane,isSpectrometerPlane)) besttc = planeClusterItr.Ptr();
                }
                planeClusterItr.Next();
              }//end loop over the clusters in plane currentPlane
              
              //check to see if there are any clusters which need to be removed
              //from the fit between the plane of planes before and that of 
              //planeCluster with the addition of planeCluster the 
              //CheckForBadClusters returns a true if it is ok to remove bad clusters
              //and zero's besttc if there were not - so if it is true,
              //then we add besttc to the track
              if(CheckForBadClusters(besttc, track, planesBefore, direction, 
                                     trk2dmin) ){
                
                MSG("TrackSR", Msg::kDebug) << " adding tc "
                                            << besttc->GetPlane()
                                            << " " 
                                            << besttc->GetTPos()
                                            << endl;
                AddBestPlaneClusterToTrack(besttc,track);
                
                //reset hitPlanesSkippedInTrack because you just added a hit plane to the track
                hitPlanesSkippedInTrack = 0;
              }
              //no besttc on this plane
              else if(direction*currentPlane
                      >direction*track->GetCluster(track->GetLast())->GetPlane())
                ++hitPlanesSkippedInTrack;
              
            }//end if distances between planes is reasonable
              
            //the distances were not reasonable, so why waste any more 
            //time on this track by making the distance even larger?  
            //set planesBefore to 999999 to break out of the loop for 
            //this track
            else planesBefore = 999999;
            
            MSG("TrackSR", Msg::kDebug) << "planes before = " << planesBefore 
                                        << " track->GetLast() = " 
                                        << track->GetLast()+1
                                        << endl;
            ++planesBefore;
          }//end loop to look for a cluster on this plane
        }//end if currentPlane is downstream of the beginning plane for this loop  
        
        //if you have a stupid number for planes before, stop the for loop and move on to the next track
        if(planesBefore >= 999999) break;
      }//end loop over planes hit in slice
    }//end if track made in this iteration

    ++trackCtr;
    
  }//end loop over tracks

  trackItr.ResetFirst();
  // check each track for outlier hit at two track ends, and remove if found
  while( (track = trackItr()) ){
    if(track->GetLast()>2){
      Int_t gap1=TMath::Abs(track->GetCluster(0)->GetPlane()-track->GetCluster(1)->GetPlane())/TMath::Abs(track->GetCluster(2)->GetPlane()-track->GetCluster(1)->GetPlane());
      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());
      Int_t gapfactor=2;
      if(gap1>fParmTrk2DNSkipHit*gapfactor)track->RemoveAt(0);
      if(gap2>fParmTrk2DNSkipHit*gapfactor)track->RemoveAt(track->GetLast());
      track->Compress();
    }   
  }
  trackItr.ResetFirst();
  return hitPlanesSkippedInTrack;
}

//_______________________________________________________________
// method to determine whether this u/v position should is active

Bool_t AlgTrackSRList::PlaneIsActive(Int_t plane, Float_t u, Float_t v, Float_t projErr)
{
  UgliGeomHandle ugh(fVldc); 
  UgliScintPlnHandle planeid;  
  PlexPlaneId plexPlane(fDetector,plane, 0); 
  if(!plexPlane.IsValid()) return false; 
  if(plexPlane.GetPlaneCoverage() == PlaneCoverage::kNoActive) return false;
  if(projErr<0.3)projErr=0.3;
  float x = 0.707*(u-v);
  float y = 0.707*(u+v);
  float dist,xedge,yedge;
  fPL.DistanceToEdge(x, y,
                    plexPlane.GetPlaneView(),
                    plexPlane.GetPlaneCoverage(),
                    dist, xedge, yedge);
  Bool_t isInside = fPL.IsInside( x, y,
                    plexPlane.GetPlaneView(),
                   plexPlane.GetPlaneCoverage());
  isInside &= dist>projErr;
  return isInside;
       
}

//----------------------------------------------------------------------------
//method to find the number of planes between currentPlane and prevPlane 
//which are skipped and active. dont count planes that are either not active, 
//the track goes through the coil hole in the plane, or the plane contains a 
//wide cluster
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)
{ 
  UgliGeomHandle ugh(fVldc);

  PlaneView::PlaneView_t view0 = PlaneView::kU;
  PlaneView::PlaneView_t view1 = PlaneView::kV;

  Int_t thisPlane = prevPlane+direction;
  Int_t numSkipped = 0;

  //projected transverse position at the currentPlane
  //  Float_t tPosProj = 999.;
  Float_t tPosProj = prevTPos;
  UgliScintPlnHandle planeid;
  PlexPlaneId plexPlane(fDetector,prevPlane, 0);
  PlexPlaneId plexPlaneEnd(fDetector,currentPlane, 0);

  plexPlane = plexPlane.GetAdjoinScint(direction);

  while(plexPlane != plexPlaneEnd && numSkipped<=fParmTrk2DNSkip){
    
    thisPlane = plexPlane.GetPlane();

    MSG("TrackSR",Msg::kDebug) << "looking at plane " << thisPlane 
                              << " view = " << plexPlane.GetPlaneView() 
                              << endl;

    // !! Consider as "spectrometer planes" the ones that are 
    // in the partially instrumented region of the detector, in order 
    // to better accomodate track finding there

    if(plexPlane.IsValid()  
       && ( (plexPlane.GetPlaneCoverage() != PlaneCoverage::kNoActive 
             &&  fDetector==Detector::kFar) 
            || (plexPlane.GetPlaneCoverage()==PlaneCoverage::kNearPartial 
                && ((plexPlane.GetPlaneView()==view0 && (tPosProj > 0. && tPosProj<2.4)) 
                    || (plexPlane.GetPlaneView()==view1 && (tPosProj  < 0. && tPosProj>-2.4)))) 
            || (plexPlane.GetPlaneCoverage()==PlaneCoverage::kNearFull)) 
       && plexPlane.GetPlaneView() == view)     
      {           
        
        planeid = ugh.GetScintPlnHandle(plexPlane);
        
        tPosProj = prevTPos;
        //      tPosProj = 999.;
        
        if(planeid.IsValid()) 
          tPosProj = prevTPos + (planeid.GetZ0()-prevZPos)*trackSlope;
        
        MSG("TrackSR",Msg::kDebug) << "projected position " << tPosProj 
                                   << " iswide plane = " 
                                   << fMapIsWide[thisPlane] 
                                   << "  dplaneoff = " << numSkipped 
                                   << " tpos = " << prevTPos 
                                   << "  slope = " << trackSlope 
                                   << endl;
      
        //increment numSkipped if any of the above reasons are true
        if(planeid.IsValid() &&  fMapIsWide[thisPlane]==0 
           && ((fDetector==Detector::kFar
                && TMath::Abs(tPosProj)>0.3) 
               || (fDetector==Detector::kNear 
                   && TMath::Abs(tPosProj)>0.8))
           ) 
          ++numSkipped;
        
      }//end if plexplane is valid
    
    plexPlane = plexPlane.GetAdjoinScint(direction);
  }//end while looking for missing planes

  return numSkipped;
}
//----------------------------------------------------------------------------
//method to find the number of planes between currentPlane and prevPlane 
//which are skipped and active. dont count planes that are either not active, 
//the track goes through the coil hole in the plane, or the plane contains a 
//wide cluster
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)
{ 
  UgliGeomHandle ugh(fVldc);

  Int_t thisPlane = prevPlane+direction;
  Int_t numSkipped = 0;

  //projected transverse position at the currentPlane
  //  Float_t tPosProj = 999.;

  Float_t uPosProj = prevUPos;
  Float_t vPosProj = prevVPos;
  UgliScintPlnHandle planeid;
  PlexPlaneId plexPlane(fDetector,prevPlane, 0);
  PlexPlaneId plexPlaneEnd(fDetector,currentPlane, 0);
  plexPlane = plexPlane.GetAdjoinScint(direction);

  // break out of loop when count exceeds 5 times value of nskip
  // parameter (max # skipped planes allowed in spectrometer)

  while(plexPlane != plexPlaneEnd && numSkipped<=fParmTrk2DNSkip*5){
    
    thisPlane = plexPlane.GetPlane();
    //increment numSkipped if any of the above reasons are true
    if(plexPlane.IsValid() &&  fMapIsWide[thisPlane]==0 
       && plexPlane.GetPlaneView()==view) {
      MSG("TrackSR",Msg::kDebug) << "looking at plane " << thisPlane 
                                 << " view = " << plexPlane.GetPlaneView() 
                                 << " coverage " << plexPlane.GetPlaneCoverage() << endl;
      
      // !! Consider as "spectrometer planes" the ones that are 
      // in the partially instrumented region of the detector, in order 
      // to better accomodate track finding there
      
        planeid = ugh.GetScintPlnHandle(plexPlane);
        uPosProj = prevUPos;
        vPosProj = prevVPos;
        if(planeid.IsValid()){ 
          uPosProj = prevUPos + (planeid.GetZ0()-prevZPos)*trackSlopeU;
          vPosProj = prevVPos + (planeid.GetZ0()-prevZPos)*trackSlopeV;
          
            MSG("TrackSR",Msg::kDebug) << " projected U position " << uPosProj 
                             << " projected V position " << vPosProj 
                             << " prev U position " << prevUPos
                             << " prev V position " << prevVPos
                             << " iswide plane = " 
                             << fMapIsWide[thisPlane] 
                             << "  dplaneoff = " << numSkipped << endl 
                             << "  slope = " << trackSlopeU << "/" 
                             << trackSlopeV 
                             << endl;
          
          Float_t posErr = TMath::Abs(planeid.GetZ0()-prevZPos)*fParmExtrapError;
          if(posErr>0.2)posErr=0.2;
          if( PlaneIsActive(thisPlane,uPosProj,vPosProj, posErr)){
            ++numSkipped;
            MSG("TrackSR",Msg::kDebug) << "plane is active " << numSkipped << endl;
          }
        } // end if planeid is valid 
      }//end if plexplane is valid
    
    plexPlane = plexPlane.GetAdjoinScint(direction);
  }//end while looking for missing planes
  
  return numSkipped;
}

//----------------------------------------------------------------------------
//method to find the number of hit planes between currentPlane and prevPlane 
Int_t AlgTrackSRList::FindNumSkippedPlanes(Int_t currentPlane, Int_t prevPlane,
                                           std::vector<Int_t>& hitPlanes, 
                                           Int_t direction)
{ 
  Int_t numSkipped = 0;

  for(Int_t i = 0; i < static_cast<Int_t>(hitPlanes.size()); ++i){
    
    if(direction*hitPlanes[i]>direction*prevPlane
       && direction*hitPlanes[i]<direction*currentPlane) ++numSkipped;
  }

  return numSkipped;
}
//---------------------------------------------------------------------
//method to identify if this is the best cluster in the plane
//the residual parameters are tested against those in the residParams
//array and if they are better, we reset the values in the array and 
//return a true
Bool_t AlgTrackSRList::IsBestClusterInPlane(TrackClusterSR *planeCluster,
                                            TrackClusterSR *prevCluster, 
                                            Int_t numSkippedHit, Int_t nHit,
                                            Int_t nBef,
                                            Double_t trackSlope, Track2DSR *track, 
                                            Double_t *residParams, 
                                            Int_t dplane, 
                                            Bool_t isSpectrometerPlane)
{
  
  bool bestCluster = false;
  bool useph = (fDetector == Detector::kFar && fParmIsCosmic==1);

 // various constants used in this method
  Double_t sigmams = 0.014;
  Double_t radLenPerPlane = 1.4;
  Double_t dPPerPlane=40*Munits::MeV;
  Double_t nomBField=0.8; // Telsa

  //get a pointer to a cluster
  TrackClusterSR *tctmp = 0;

  //get variables to describe the differences from the prevCluster to the 
  //one in the current plane in z, time, transverse position,etc
  Double_t dzpos = TMath::Abs(planeCluster->GetZPos()-prevCluster->GetZPos());
  Double_t dtime = planeCluster->GetBegTime()-track->GetT0()-
    (planeCluster->GetZPos()-track->GetZ0())/Mphysical::c_light;
   
  Double_t resid=0.;
  Double_t maxresid=99999.;
  Double_t dangle = 0.;
  Double_t maxdangle = 999999.;

  //declare the size of the window in # of planes over which you want to 
  //find the slope of the track
  // init maximum allowable residual fixed sized (related to strip width) plus
  // contribution scaling with distance extrapolated
  Int_t mhit;
  if(nHit<fParmTrk2DNHough0){
    mhit = nHit;    
    maxresid = fParmTrk2DLinA0+fParmTrk2DLinB0*TMath::Abs(dzpos);
  }
  else{
    mhit = min(nHit,fParmTrk2DNHough);
    if( fDetector == Detector::kNear && prevCluster->GetPlane()>121 && mhit>2)mhit--;
    maxresid = TMath::Abs(planeCluster->GetMaxTPos()-planeCluster->GetMinTPos())*0.288;
  }
  maxresid *=maxresid;
  if(nHit>1){     

    //nexitplane tells you how many planes you have left between the 
    //plane of the last cluster added to the track and the most downstream hit in the
    //slice. This is a (crude) estimate of the muon energy at that point,  and
    // allows an estimate of bending and MCS 
    Int_t nexitplane = TMath::Abs(prevCluster->GetPlane()-fSliceEndPlane)+1;
    if(fSliceNUPlanes+fSliceNVPlanes<fParmTrkNPlaneGoodDir && fParmIsCosmic==1){
      // this is a short track, so the timing is not so good
      // take smaller momentum from either end
      nexitplane = min(nexitplane,TMath::Abs(prevCluster->GetPlane()-fSliceVtxPlane)+1);
    }
       
    Double_t *xfit = new Double_t[mhit];
    Double_t *yfit = new Double_t[mhit];
    Double_t *wfit = new Double_t[mhit];
    Double_t parm[2];
    Double_t eparm[2] = {0.,0.};
    Double_t wtmp = 0.;
 
    //fit a straight line over the clusters in a window of size
    //mhit to see how well the current cluster matches up to the
    //fit of previously added clusters
    for (Int_t ifit=0; ifit<mhit; ++ifit) {
      if(isSpectrometerPlane){
        tctmp = track->GetCluster(ifit+nBef);
      }
      else{
        tctmp = track->GetCluster(nHit-mhit+ifit);
      }
      xfit[ifit] = tctmp->GetZPos()-planeCluster->GetZPos();
      yfit[ifit] = tctmp->GetTPos();
      wtmp = 0.288*(tctmp->GetMaxTPos()-tctmp->GetMinTPos());
      wfit[ifit]=1.0/(wtmp*wtmp);
    }

    for(Int_t ifit=0;ifit<mhit;ifit++){
      MSG("TrackSR", Msg::kDebug)  << " fit input " << xfit[ifit]  << " " << yfit[ifit] << " " << wfit[ifit] << endl; 
    }
    
    LinearFit::Weighted(mhit,xfit,yfit,wfit,parm,eparm);
    
    delete [] xfit;
    delete [] yfit;
    delete [] wfit;

    //get the residual for this plane to the fit over the previous planes
    MSG("TrackSR", Msg::kDebug) << "fit slope = " << parm[1] <<  "/ " << eparm[1] 
                                << " fit intercept = " 
                                << parm[0] << " / " << eparm[0] << endl;

    //find the change in the angle from the fit slope 
    Double_t planeClusterSlope = (planeCluster->GetTPos()-prevCluster->GetTPos())
                                  /(planeCluster->GetZPos()-prevCluster->GetZPos());

    dangle = TMath::Abs(parm[1]-planeClusterSlope);

    resid = TMath::Abs(parm[0]-planeCluster->GetTPos());

    // determine max allowable residuals
    //get the range-based  momentum assuming a loss of 40 MeV per plane
    Double_t dzds = 1.0/TMath::Sqrt(1.+trackSlope*trackSlope);
    Double_t pathlength = (Double_t)(nexitplane)/dzds;
    Double_t tctp = pathlength*dPPerPlane;
    Double_t dbfangle =  0.3*dzpos/dzds*nomBField/tctp;
    Double_t mcsangle =  sigmams*TMath::Sqrt(radLenPerPlane*TMath::Abs(dplane)/dzds)/tctp;     

    maxdangle = planeCluster->GetTPosError()*planeCluster->GetTPosError()/(dzpos*dzpos) +
      dbfangle*dbfangle + mcsangle*mcsangle + eparm[1]*eparm[1] ; 
    maxresid += dzpos*dzpos*(dbfangle*dbfangle + mcsangle*mcsangle) +  eparm[0]*eparm[0]; 
 
    maxresid  =  fParmTrk2DNSigmaA*TMath::Sqrt(maxresid);
    maxdangle =  fParmTrk2DDirCosNSigma*TMath::Sqrt(maxdangle);

    // make requirements more stringent if skipping already found clusters
    if(nBef>0){
      maxresid /=2.;
      maxdangle /=2.;
    }
 
    MSG("TrackSR",Msg::kDebug) << " max resids:  dzpos " << dzpos 
                               << " bfield " << dbfangle 
                               << " mcs " << mcsangle 
                               << " maxdangle " << maxdangle 
                               << " bfield resid " << dbfangle*dzpos 
                               << " mcs resid " << mcsangle*dzpos << endl;    
  }//end if at least 2 hits
  
  else if(track->GetHoughExist()){
    //end not enough hits, try getting residuals using the hough transform values 
    tctmp = track->GetCluster(0);
    resid = min(TMath::Abs((tctmp->GetTPos()+(planeCluster->GetZPos()-tctmp->GetZPos())*track->GetHoughSlope())-planeCluster->GetTPos()),
                TMath::Abs(tctmp->GetTPos()-planeCluster->GetTPos()));
   
    maxresid  =  fParmTrk2DNSigmaA*TMath::Sqrt(maxresid +
                                               (tctmp->GetMaxTPos()-tctmp->GetMinTPos())*
                                               (tctmp->GetMaxTPos()-tctmp->GetMinTPos())*0.083);  
  }

  //get the weight residual based on an empirical parameterization done by
  //R. Lee
  Double_t weightedresid = (0.3+exp(-0.2*planeCluster->GetCharge()))*resid;
  maxresid = min(0.4,maxresid);  // don't allow maxresid to exceed 0.4 m
  
  MSG("TrackSR",Msg::kDebug) << "plane " << planeCluster->GetPlane() << "numskiphit " <<
    numSkippedHit << "/" << fParmTrk2DNSkipHit << endl <<
    " nhit " << nHit << "/" << fParmMinSingleHit << " nstrip " << planeCluster->GetNStrip() << "/" << fParmSingleHitDef <<  endl <<
    " dtime " << dtime*1.e9 << "/" << fParmHitNTime*1.e9 << "/" << fParmHitTime*1.e9 <<  endl << 
    " resid " << resid << "/" << maxresid  << " dangle " << dangle << "/" << maxdangle << endl <<
    " useph " << useph << "/" << residParams[0]  << "/" << TMath::Abs(weightedresid) << "/" << residParams[1] << endl; 


  // criteria for this cluster being the best cluster on this plane to add to track
  if(numSkippedHit<=fParmTrk2DNSkipHit 
     && (nHit>fParmMinSingleHit || planeCluster->GetNStrip()<=fParmSingleHitDef) 
     && TMath::Abs(resid)<maxresid && dangle<maxdangle
     && dtime>fParmHitNTime && dtime<fParmHitTime 
     && ((!useph && TMath::Abs(resid)<residParams[0]) 
         || (useph && TMath::Abs(weightedresid)<residParams[1]))){
    MSG("TrackSR",Msg::kDebug) << "  accepted with resid = " <<  resid
                               << " weighted resid = "<< weightedresid << endl;
    // add hit to this track
    bestCluster = true;
    residParams[0] = TMath::Abs(resid);
    residParams[1] = TMath::Abs(weightedresid);
  }//end if this cluster works
  
  return bestCluster;
}


//------------------------------------------------------------------------------
//this method looks at the track to see if there are any clusters between the
//last good cluster in the track and the current plane with an identified cluster
//to add to the track which should be removed from the track.  ie you have a 
//situation where a cluster was added to the track in a previous iteration, but 
//the information from the current plane makes it apparent that the tracker 
//screwed up and shouldnt have added that cluster.
Bool_t AlgTrackSRList::CheckForBadClusters(TrackClusterSR *besttc, 
                                           Track2DSR *track,
                                           Int_t &planesBefore, Int_t direction, 
                                           Int_t trk2dmin)
{
  
  Bool_t removetc(kTRUE);

  if( !besttc) return false;
  if(planesBefore<=0 ) return true;
  
   MSG("TrackSR", Msg::kDebug) << " ----------checking besttc on plane------------ " << besttc->GetPlane()
                              << endl;

   //get the number of clusters (ie hit planes) in the track.
   //add one to track->GetLast() because that returns a number starting
   //from 0 - its an array count after all, and add one to planesBefore
   //to account for besttc just found but not yet added
   Int_t nhit = track->GetLast()+1-planesBefore+1;
   
   //assume alternate views when finding dplane - the number of planes
   //between the first plane in the track and the plane of the current cluster.
   //not the best way since it has problems in tracks that jump the super module
   //gap, but good enough for now
   Int_t dplane = direction*(besttc->GetPlane()-track->GetBegPlane())/2+1;
   
   MSG("TrackSR", Msg::kDebug) << "nhit = " << nhit << " dplane = " << dplane
                               << "/" << fParmTrk2DHitFraction << " trk2dmin = "
                               << trk2dmin << " 2DNHough = " << fParmTrk2DNHough
                               << endl;
   
   //if the # of hits in the track is less than the hit fraction * dplane
   //you cant remove any of the hits - you need them all
   //if the # of hits is less than the minimum you cant get rid of any 
   if ((1.*nhit)<fParmTrk2DHitFraction*dplane || nhit<trk2dmin) {
     removetc = kFALSE;
     planesBefore = 9999999;
   }
   
   //if you only have enough hits to do hough tracks, see how well
   //the current hit fits in with the rest
   else if(nhit<=fParmTrk2DNHough){
     Double_t *xfit = new Double_t[nhit];
     Double_t *yfit = new Double_t[nhit];
     Double_t parm[2];
     TrackClusterSR *tctmp = 0;
     for (Int_t ifit=0; ifit<nhit; ++ifit){
       tctmp = track->GetCluster(ifit);
       if (ifit==nhit-1) tctmp = besttc;
       xfit[ifit] = tctmp->GetZPos();
       yfit[ifit] = tctmp->GetTPos();
     }
     LinearFit::Unweighted(nhit,xfit,yfit,parm);
     Double_t maxresid = 0.;
     for(Int_t ifit=0; ifit<nhit; ++ifit){
       maxresid = max(maxresid,
                      TMath::Abs(parm[0]+parm[1]*xfit[ifit]-yfit[ifit]));
     }
     delete [] xfit;
     delete [] yfit;
     MSG("TrackSR", Msg::kDebug)<< "maxresid = " << maxresid << " / " 
                                << fParmTrk2DMaxResid<< endl;
     if (maxresid>fParmTrk2DMaxResid) {
       removetc = kFALSE;
       planesBefore = 9999999;
     }
   }//end if not enough hits for anything but a hough track
  if(removetc){
    // solution found, need to remove bad track cluster(s) between
    //the plane of the besttc and the plane represented by
    //planesBefore
    for(int iremove=0; iremove<planesBefore; ++iremove) {
           MSG("TrackSR", Msg::kDebug)<< "removing cluster on plane  " << track->GetCluster(track->GetLast())->GetPlane() << endl;
    track->RemoveAt(track->GetLast());
      track->Compress();
    }
  } 
  else{
    //adding this cluster will somehow make the fit worse, so dont add it
    besttc = 0;
  }
  return removetc;
}

//------------------------------------------------------------------------
//method which actually adds the best cluster in a plane to the track
void AlgTrackSRList::AddBestPlaneClusterToTrack(TrackClusterSR *besttc,
                                                Track2DSR *track)
{
  
  Int_t iclosest=-1;
  Float_t dzclosest=1e6;
  Int_t ilast = track->GetLast();
  TrackClusterSR *tc0 = track->GetCluster(ilast);
  for (Int_t i=0; i<track->GetLast(); i++){
    tc0 = track->GetCluster(i);
    if(TMath::Abs(tc0->GetZPos()-besttc->GetZPos())<dzclosest){
      iclosest=i;
      dzclosest=TMath::Abs(tc0->GetZPos()-besttc->GetZPos());
    }
  }
  Double_t slope0=0;
  if(iclosest>=0){
    slope0 = track->GetForwardSlope(iclosest);
  }
  Double_t slope = (besttc->GetTPos()-tc0->GetTPos()); 
  slope /= (besttc->GetZPos()-tc0->GetZPos());
  
  //weight the slope of the this cluster using the previous clusters in the track
  if (slope0!=0) {
    slope *= (1.-fParmTrk2DAlpha);
    slope += fParmTrk2DAlpha*slope0;
  }
  track->Add(besttc,slope);
  return;
}
     
//-----------------------------------------------------------------------------
//method to id and remove bad tracks based on the hit fraction, track length, etc
  Int_t AlgTrackSRList::IdentifyBadTracks(Track2DSRItr &trackItr, Int_t iterate,
                                        Int_t trk2dmin)
{
  //get to the first track in the list
  trackItr.ResetFirst();
  Track2DSR *track = 0;
  TrackClusterSR *tctmp = 0;
  Double_t hitfrac = 0.;
  Int_t ilast = 0;
  Int_t nhit = 0;
  Int_t removedTracks = 0; // Double_t dplane = 0.;
  while( (track = trackItr()) ){ //loop over the tracks

    //if the track was created in this iteration and isnt already bad
    if(track->fIterate==iterate && !track->IsBad() ){
      hitfrac = 0.;
      Bool_t hitRemoved=true;
      while(hitRemoved && track->GetLast()>0 ){
        ilast = track->GetLast();
        nhit = ilast+1;
        hitRemoved=false;
        Float_t aveSpacing= 0.5*TMath::Abs(track->GetEndPlane()-
                                   track->GetBegPlane())/ilast;
        Float_t lastSpacing=0.5*TMath::Abs(track->GetCluster(ilast)->GetPlane()-
                                             track->GetCluster(ilast-1)->GetPlane());
        hitfrac=aveSpacing/lastSpacing;

        
        //remove the last clusters from tracks with too few hits for the number 
        //of planes crossed until the hit occupancy is acceptable
        if(hitfrac<fParmTrk2DHitFraction && lastSpacing>5){
          track->RemoveAt(ilast);
          hitRemoved=true;
          MSG("TrackSR",Msg::kDebug) << "removing end hit: hit fraction " 
                                     << hitfrac << "/" << fParmTrk2DHitFraction 
                                     <<endl;
        }

        
      }//end loop while not enough hit occupancy
        
      //remove short tracks
      if(track->GetLast()+1<trk2dmin){
        ++removedTracks;
        track->SetIsBad(true);
        MSG("TrackSR",Msg::kDebug) << "adding badtrack: short track " 
                                   << track->GetLast()+1 << "/" << trk2dmin 
                                   << " beg plane = " << track->GetBegPlane()
                                   << endl;
      }//end if short track


      //check linearity of short tracks
      nhit = track->GetLast()+1;
      if(nhit<=fParmTrk2DNHough && nhit>2){
        Double_t *xfit = new Double_t[nhit];
        Double_t *yfit = new Double_t[nhit];
        Double_t parm[2];
        for (Int_t ifit=0; ifit<nhit; ++ifit) {
          tctmp = track->GetCluster(ifit);
          xfit[ifit] = tctmp->GetZPos();
          yfit[ifit] = tctmp->GetTPos();
        }
        LinearFit::Unweighted(nhit,xfit,yfit,parm);
        Double_t maxresid = 0.;
        for (Int_t ifit=0; ifit<nhit; ++ifit) {
          maxresid = max(maxresid, TMath::Abs(parm[0]+parm[1]*xfit[ifit]-yfit[ifit]));
        }
        delete [] xfit;
        delete [] yfit;
        if(maxresid>fParmTrk2DMaxResid*max(fSliceDzDs,
                                           TMath::Sqrt(1.+parm[1]*parm[1]))){
          ++removedTracks;
          track->SetIsBad(true);
          MSG("TrackSR",Msg::kDebug) << "adding badtrack: linearity " 
                                     << maxresid << "/" 
                                     << fParmTrk2DMaxResid <<endl;
        }

      }//end if short track

    }//end if track created in current iteration

    MSG("TrackSR", Msg::kDebug) << "track with beg plane " << track->GetBegPlane()
                               << " is bad " << (Int_t)track->IsBad() << endl;
  }//end loop over tracks
 
  return removedTracks;
}

//----------------------------------------------------------------------
//this method marks all clusters used in the good tracks from the
//current iteration as non-valid so that they are not used in new tracks
void AlgTrackSRList::MarkUsedClusters(Track2DSRItr &trackItr, 
                                      TrackClusterSRItr &clusterItr, 
                                      Int_t iterate, Int_t nmaxtrack,
                                      Double_t houghSlope, Double_t houghIntercept)
{
  Track2DSR *track = 0;
  TrackClusterSR *tc = 0;
  TrackClusterSR *tc1 = 0;
  TrackClusterSR *tcbeg = 0;
  TrackClusterSR *tcend = 0;

  Double_t residbeg =0.;
  Double_t residend =0.;

  trackItr.ResetFirst();
  clusterItr.ResetFirst();

  //loop over the tracks
  while( (track = trackItr()) ){
    
    //if this track was created in this iteration
    if(track->fIterate==iterate){

      //if it was a good track
      if( !track->IsBad() ){
        MSG("TrackSR",Msg::kDebug) << "removing hits from track\n";

        //loop over the clusters in this track
        for(Int_t itc=0; itc<=track->GetLast(); ++itc){
          tc = track->GetCluster(itc);
          MSG("TrackSR",Msg::kDebug) << "  TC " << tc->GetPlane() << " " 
                                       << tc->GetTPos() <<endl;

          //loop over all clusters in the clusterItr and mark all the ones
          //used in this track as non-valid
          while( (tc1 = clusterItr()) ){

            //if tc1 is valid and the same cluster as tc, mark it as non-valid
            //because it has been used
            if(tc1->IsValid() && tc1->GetPlane()==tc->GetPlane() 
               && TMath::Abs(tc1->GetTPos()-tc->GetTPos())<1.e-3
               && TMath::Abs(1.e9*(tc1->GetBegTime()-tc->GetBegTime()))<1.e-3){
              
              tc1->SetIsValid(false);

              MSG("TrackSR",Msg::kDebug) << "    REMOVING\n";
            }//end if the same cluster and valid

          }//end loop over clusters

          clusterItr.ResetFirst();
        }//end loop over clusters in this track
      }//end if good track
      else{ //now do the bad tracks

        if(nmaxtrack>0){
          tcbeg = track->GetCluster(0);
          tcend = track->GetCluster(track->GetLast());
     
          residbeg = TMath::Abs(houghIntercept
                                +houghSlope*tcbeg->GetZPos()-tcbeg->GetTPos());
          residend = TMath::Abs(houghIntercept
                                +houghSlope*tcend->GetZPos()-tcend->GetTPos());
          MSG("TrackSR",Msg::kDebug) << "bad tracks, hough track exists: tcbeg " 
                                     << tcbeg->GetPlane() << " " 
                                     << tcbeg->GetTPos() 
                                     << " " << residbeg << "   tcend " 
                                     << tcend->GetPlane() << " " 
                                     << tcend->GetTPos() 
                                     << " " << residend << endl;

          //remove the cluster with the larger residual
          if(residbeg<residend){
            MSG("TrackSR",Msg::kDebug) << "  removing end cluster\n";
            tc = tcend;

            //loop over all clusters in the clusterItr and mark all the ones
            //used in this track as non-valid
            while( (tc1 = clusterItr()) ){
            
              if (tc1->IsValid() && tc1->GetPlane()==tc->GetPlane() &&
                  TMath::Abs(tc1->GetTPos()-tc->GetTPos())<1.e-3 &&
                  TMath::Abs(1.e9*(tc1->GetBegTime()-tc->GetBegTime()))<1.e-3) {
                tc1->SetIsValid(false);
                MSG("TrackSR",Msg::kDebug) << "  REMOVING HIT\n";
              }//end if the same cluster
            }//end loop over clusters
            clusterItr.ResetFirst();
            
          } //end if removing end cluster
          else{ //removing begining cluster
            MSG("TrackSR",Msg::kDebug) << "  removing beg cluster\n";
            tc = tcbeg;
          }
        }//end if at least one track in this view
        else{
          tc = track->GetCluster(0);
          MSG("TrackSR",Msg::kDebug) << "bad tracks, no hough track: tcbeg " 
                                       << tc->GetPlane() << " " << tc->GetTPos() 
                                       << " " 
                                       << tc->GetBegTime()*1.e9 << "\n";
        }
        
        //loop over all clusters in the clusterItr and mark all the ones
        //used in this track as non-valid
        while( (tc1 = clusterItr()) ){
          if (tc1->IsValid() && tc1->GetPlane()==tc->GetPlane() &&
              TMath::Abs(tc1->GetTPos()-tc->GetTPos())<1.e-3 &&
              TMath::Abs(1.e9*(tc1->GetBegTime()-tc->GetBegTime()))<1.e-3) {
            tc1->SetIsValid(false);
            MSG("TrackSR",Msg::kDebug) << "  REMOVING HIT\n";
          }//end if the same cluster
        }//end loop over clusters
        clusterItr.ResetFirst();

      }//end else bad track
    }//end if track created in this iteration
  }//end loop over tracks
  
  return;
}

//----------------------------------------------------------------
//this method removes track subsets to get rid of redundant tracks
Int_t AlgTrackSRList::RemoveTrackSubsets(Track2DSRItr &trackItr1, 
                                         Track2DSRItr &trackItr2, 
                                         TObjArray *trackList)//, Int_t direction)
{
  Track2DSR *track1 = 0;
  Track2DSR *track2 = 0;

  trackItr1.ResetFirst();
  trackItr2.ResetFirst();

  Int_t removeCtr = 0;

  //counters to keep track of which track we are on in 
  //each set - dont want to remove a track because we happen to
  //be on the same one in the different sets

  Int_t itrk1 = 0;
  Int_t itrk2 = 0;
 
  Double_t chi21 = 0.;
  Double_t chi22 = 0.;
  
  Int_t nmatch = 0;

  TrackClusterSR *tc1 = 0;
  TrackClusterSR *tc2 = 0;

  Bool_t found = true;
  Bool_t match = false;
  Double_t matchfraction = 0.;

  //loop over tracks in both iterators to see if any match 
  while( (track1 = trackItr1()) ){
    
    //if there are any clusters in this track and it isnt a bad track
    if (track1->GetLast()>=0 && !track1->IsBad() ) {

       MSG("TrackSR", Msg::kDebug) << "track 1 " << "begPlane = " 
                                   << track1->GetBegPlane()
                                   << " end plane = " << track1->GetEndPlane() 
                                   << " itrk1  " << itrk1 << " last is  " 
                                   << track1->GetLast()
                                   << endl;


      itrk2 = 0;
      while( (track2 = trackItr2()) ){
    
        MSG("TrackSR", Msg::kDebug) << "track 2 " << "begPlane = " 
       
                                   << track2->GetBegPlane()
                                   << " end plane = " << track2->GetEndPlane() 
                                   << " itrk2 " << itrk2 << " last is  " << track2->GetLast() 
                                   << endl;
      
        if(track2->GetLast()>=0 && !track2->IsBad()){
          //if not on the same track but track2 has fewer clusters than
          //track1
          //if track 2 begins before track1 ends
          if( itrk1!=itrk2 
              && track2->GetLast()<=track1->GetLast() ){
                    
            MSG("TrackSR", Msg::kDebug) << "track2 could be a subset of track 1" 
                                        << endl; 
                
            // check if track2 is subset of track1
            found = true;
            nmatch=0;
            
            //loop over the clusters in the 2 tracks to see if any are the same
            for(Int_t itc2=0; itc2<=track2->GetLast(); ++itc2){
              tc2 = track2->GetCluster(itc2);
              match = false;
              for(Int_t itc1=0; itc1<=track1->GetLast() && !match; ++itc1){
                tc1 = track1->GetCluster(itc1);
                if(tc1->GetPlane()==tc2->GetPlane() 
                   && tc1->GetMinStrip()==tc2->GetMinStrip() ){
                  match = true;
                  MSG("TrackSR" , Msg::kDebug) << "matching cluster found" << endl; 
                }
              }//end loop over track1 clusters

              //if none are similar, no subset found
              if(!match){
                MSG("TrackSR" , Msg::kDebug) << "matching cluster not found" 
                                             << endl;
                found = false;
              }
              else ++nmatch;
              
            }//end loop over track2 clusters
                  
            //if found is still set then all the clusters in track2 were 
            //in track 1, so set track2 as bad
            if(found){
              track2->SetIsBad(true);
              MSG("TrackSR", Msg::kDebug) << "all track2 clusters found in track1 "
                                          << "set track " << itrk2 << " bad" 
                                          << endl;
            }
            else{
              matchfraction = (1.*nmatch)/(1.*track2->GetLast()+1.);
             MSG("TrackSR" , Msg::kDebug) << "matchfraction = " << matchfraction 
                                          << " twin match fraction = " 
                                          << fParmTrk2DTwinMatchFraction 
                                          << " nmatch = " << nmatch 
                                          << " subset n hits = " << fParmTrk2DSubsetNHit 
                                          << endl; 

              if(track2->GetLast()+1-nmatch<fParmTrk2DSubsetNHit && nmatch>0){
                track2->SetIsBad(true);
                MSG("TrackSR", Msg::kDebug) << "setting track2 bad" << endl;
                
              }
              else if(matchfraction>fParmTrk2DTwinMatchFraction){
                chi21 = track1->GetChi2()/(1.*track1->GetLast()+1.);
                chi22 = track2->GetChi2()/(1.*track2->GetLast()+1.);
                MSG("TrackSR", Msg::kDebug) << "chi21 = " << chi21
            
                                           << " chi22 = " << chi22 << endl;
                if(chi21>chi22){
                  track1->SetIsBad(true);
                  MSG("TrackSR", Msg::kDebug) << "setting track1 bad, chi2" << endl;
                  
                }
                else{
                  track2->SetIsBad(true);
                  MSG("TrackSR", Msg::kDebug) << "setting track2 bad, chi2" << endl;
                 
                }
              }//end if matchfraction makes it look like a twin
            }//end else for no matches found
          }//end if not the same track
        }//end if track2 is ok to use
        ++itrk2;
      }//end loop over trackItr2
      trackItr2.ResetFirst();
    }//end if track1 is ok to use
    ++itrk1;
  }//end loop over trackItr1


  //remove the bad tracks from the list
  for (Int_t itrk=0; itrk<=trackList->GetLast(); itrk++) {
    track1 = dynamic_cast<Track2DSR*>(trackList->At(itrk));
    if(track1->IsBad() || track1->fIterate<0){

      MSG("TrackSR" , Msg::kDebug) << "removing bad track " << itrk 
  
                                  << " with beg plane = "
                                  << track1->GetBegPlane() 
                                  << " is bad " << (Int_t)track1->IsBad()
                                  << " iteration " << track1->fIterate << endl;

      trackList->RemoveAt(itrk);
      delete track1;
      ++removeCtr;
    }
  }
  trackList->Compress();

  return removeCtr;
}
// ---------------------------------------------------------------------
// merge tracks which are 'head to tail', and consistent
void AlgTrackSRList::MergeTracks(TObjArray * tracklist, Int_t direction)

{
  MSG("TrackSR",Msg::kDebug) << " in MergeTracks" << endl;
  Track2DSRItr trackItr1(tracklist);  
  Track2DSRItr trackItr2(tracklist);
  Track2DSR *track1 = 0;
  Track2DSR *track2 = 0;
  for (Int_t nremove1=0;nremove1<5;nremove1++){ 
    for (Int_t nremove2=0;nremove2<5;nremove2++){
      Int_t fwd1,fwd2,bck1,bck2;
      Bool_t merge=false;
      Int_t itrack1=0;
      trackItr1.Reset();
      while( (track1 = trackItr1()) ){
          MSG("TrackSR",Msg::kDebug) <<  " track 1 " << track1->GetCluster(0)->GetPlane() << "/" << 
                    track1->GetCluster(track1->GetLast())->GetPlane() << endl;
        if(track1->GetLast()>fParmTrk2DNSeed && !track1->IsBad()){
          Int_t itrack2=0;
          trackItr2.Reset();
          while((track2=trackItr2())){ 
          MSG("TrackSR",Msg::kDebug) <<  " track 2 " << track2->GetCluster(0)->GetPlane() << "/" << 
                    track2->GetCluster(track2->GetLast())->GetPlane() << endl;
            if(track2->GetLast()>fParmTrk2DNSeed && !track2->IsBad()){
              if(track1->GetPlaneView()==track2->GetPlaneView()){
                merge=false;
                if(track1!=track2){
                  if(direction>0){
                    fwd1=min(nremove1,track1->GetLast()-2);
                    bck1=max(track1->GetLast()-nremove1,2);
                    fwd2=min(nremove2,track2->GetLast()-2);
                    bck2=max(track2->GetLast()-nremove2,2);
                  }
                  else{
                    fwd1=max(track1->GetLast()-nremove1,2);
                    bck1=min(nremove1,track1->GetLast()-2);
                    fwd2=max(track2->GetLast()-nremove2,2);
                    bck2=min(nremove2,track2->GetLast()-2);
                  }
                  MSG("TrackSR",Msg::kDebug) << " Tracks " << itrack1 << "/" << itrack2 << " track 1 " << track1->GetCluster(fwd1)->GetPlane() << "/" << 
                    track1->GetCluster(bck1)->GetPlane() << "    track 2 " << track2->GetCluster(fwd2)->GetPlane() << "/" << track2->GetCluster(bck2)->GetPlane() << endl;
                  if(track2->GetCluster(bck2)->GetPlane()<
                     track1->GetCluster(fwd1)->GetPlane()){  // track1 downstream of track2
                    MSG("TrackSR",Msg::kDebug) << "track 1 is downstream of track 2" << endl;
                    
                    
                    Float_t upos,vpos,uslope,vslope;
                    if(track1->GetPlaneView()==PlaneView::kU){
                      upos = track1->GetCluster(fwd1)->GetTPos();
                      vpos = fHoughVIntercept + track1->GetCluster(fwd1)->GetZPos()*fHoughVSlope;
                      uslope=track1->GetForwardSlope(fwd1);
                      vslope=fHoughVSlope;
                    }
                    else{
                      vpos = track1->GetCluster(fwd1)->GetTPos();
                      upos = fHoughUIntercept + track1->GetCluster(fwd1)->GetZPos()*fHoughUSlope;
                      vslope=track1->GetForwardSlope(fwd1);
                      uslope=fHoughUSlope;
                    }
                    Int_t dir=-1;       
                    Int_t numSkipped=FindNumSkippedPlanes(track2->GetCluster(bck2)->GetPlane(),
                                                          track1->GetCluster(fwd1)->GetPlane(),
                                                          upos,vpos,
                                                          track1->GetCluster(fwd1)->GetZPos(),
                                                          uslope,
                                                          vslope,
                                                          dir,
                                                          track1->GetPlaneView());
                    MSG("TrackSR",Msg::kDebug) << "# skipped planes = " << numSkipped << " tpos " << track1->GetCluster(fwd1)->GetTPos() << " " << track2->GetCluster(bck2)->GetTPos() << endl;
                    if((TMath::Abs(track1->GetCluster(fwd1)->GetTPos())<0.375 && TMath::Abs(track2->GetCluster(bck2)->GetTPos())<0.375  ) || numSkipped<=fParmTrk2DNSkipHit ){ 

                      Float_t midZ=0.5*(track1->GetCluster(fwd1)->GetZPos()+
                                        track2->GetCluster(bck2)->GetZPos());
                      Float_t delZ=TMath::Abs(track1->GetCluster(fwd1)->GetZPos()-track2->GetCluster(bck2)->GetZPos());
                      Float_t textrap1=track1->GetCluster(fwd1)->GetTPos()+
                        track1->GetForwardSlope(fwd1)*(midZ-track1->GetCluster(fwd1)->GetZPos());        
                      Float_t textrap2=track2->GetCluster(bck2)->GetTPos()+
                        track2->GetForwardSlope(bck2)*(midZ-track2->GetCluster(bck2)->GetZPos());
                      Float_t delTPos=TMath::Abs(textrap1-textrap2);
                      Float_t tPosErr=  TMath::Sqrt(track1->GetCluster(fwd1)->GetTPosError()*track1->GetCluster(fwd1)->GetTPosError()+
                                                    track2->GetCluster(bck2)->GetTPosError()*track2->GetCluster(bck2)->GetTPosError());   
                      Float_t delSlope= TMath::Abs(track1->GetForwardSlope(fwd1)-track2->GetForwardSlope(bck2));
                      // multiply maxResid, maxSlopeDiff by root(2) since we are extrapolating two tracks to center, and not accounting for bending, mcs, etc 
                      Float_t maxResid = 1.4*(fParmTrk2DLinA+fParmTrk2DLinB*0.5*delZ);
                      Float_t maxSlopeDiff = 1.4*(TMath::Sqrt(tPosErr*tPosErr/(delZ*delZ)+fParmTrk2DDirCosError2));
                      MSG("TrackSR",Msg::kDebug) << " delTPos " << delTPos << "/" << maxResid << " delSlope " << delSlope <<"/" << maxSlopeDiff << endl;
                      if(delTPos<maxResid   &&
                         delSlope< maxSlopeDiff ){
                        merge=true;
                        MSG("TrackSR",Msg::kDebug) << " merging tracks " << endl;
                      }
                      
                    }
                    if(merge){
                      if( nremove1 || nremove2 )  MSG("TrackSR",Msg::kDebug) << " removing hits " << nremove1 << " " << nremove2 <<  endl;
                      for(Int_t i=0;i<nremove1;i++){
                        if(direction>0){
                          track1->RemoveAt(0);
                        }
                        else{
                          track1->RemoveAt(track1->GetLast());
                        }
                        track1->Compress();
                      }
                      for(Int_t i=0;i<nremove2;i++){
                        if(direction>0){
                          track2->RemoveAt(track2->GetLast());
                        }
                        else{
                          track2->RemoveAt(0);
                        }
                        track2->Compress();
                      }
                    }
                  }
                  else if (track1->GetCluster(bck1)->GetPlane()<
                           track2->GetCluster(fwd2)->GetPlane()){ // track2 downstream of track 1
                    MSG("TrackSR",Msg::kDebug) << "track 2 is downstream of track 1" << endl;
                    
                    Int_t dir=-1;
                    Float_t upos,vpos,uslope,vslope;
                    if(track1->GetPlaneView()==PlaneView::kU){
                      upos = track2->GetCluster(fwd2)->GetTPos();
                      vpos = fHoughVIntercept + track2->GetCluster(fwd2)->GetZPos()*fHoughVSlope;
                      uslope=track2->GetForwardSlope(fwd2);
                      vslope=fHoughVSlope;
                    }
                    else{
                      vpos = track2->GetCluster(fwd2)->GetTPos();
                      upos = fHoughUIntercept + track2->GetCluster(fwd2)->GetZPos()*fHoughUSlope;
                      vslope=track2->GetForwardSlope(fwd2);
                      uslope=fHoughUSlope;
                    }
                    Int_t numSkipped=FindNumSkippedPlanes(track1->GetCluster(bck1)->GetPlane(),
                                                          track2->GetCluster(fwd2)->GetPlane(),
                                                          upos,vpos,
                                                          track2->GetCluster(fwd2)->GetZPos(),
                                                          uslope,
                                                          vslope,
                                                          dir,
                                                          track1->GetPlaneView());
                    MSG("TrackSR",Msg::kDebug) << "# skipped planes = " << numSkipped << " tpos " << track1->GetCluster(bck1)->GetTPos() << " " << track2->GetCluster(fwd2)->GetTPos() << endl; 
                    if((TMath::Abs(track1->GetCluster(bck1)->GetTPos())<0.375 && TMath::Abs(track2->GetCluster(fwd2)->GetTPos())<0.375  ) || numSkipped<=fParmTrk2DNSkipHit){             
                      Float_t midZ=0.5*(track1->GetCluster(bck1)->GetZPos()+
                                        track2->GetCluster(fwd2)->GetZPos());
                      Float_t delZ=TMath::Abs(track1->GetCluster(bck1)->GetZPos()-track2->GetCluster(fwd2)->GetZPos());
                      Float_t textrap1=track1->GetCluster(bck1)->GetTPos()+
                        track1->GetForwardSlope(bck1)*(midZ-track1->GetCluster(bck1)->GetZPos());        
                      Float_t textrap2=track2->GetCluster(fwd2)->GetTPos()+
                        track2->GetForwardSlope(fwd2)*(midZ-track2->GetCluster(fwd2)->GetZPos());
                      Float_t delTPos = TMath::Abs(textrap1-textrap2);
                      Float_t tPosErr =  TMath::Sqrt(track1->GetCluster(bck1)->GetTPosError()*track1->GetCluster(bck1)->GetTPosError()+
                                                     track2->GetCluster(fwd2)->GetTPosError()*track2->GetCluster(fwd2)->GetTPosError());          
                      Float_t delSlope = TMath::Abs(track1->GetForwardSlope(bck1)-track2->GetForwardSlope(fwd2));
                      Float_t maxResid = 1.4*(fParmTrk2DLinA+fParmTrk2DLinB*0.5*delZ);
                      Float_t maxSlopeDiff = 1.4*(TMath::Sqrt(tPosErr*tPosErr/(delZ*delZ)+fParmTrk2DDirCosError2));
                      MSG("TrackSR",Msg::kDebug) << " delTPos " << delTPos << "/" << maxResid << " delSlope " << delSlope <<"/" << maxSlopeDiff << endl;
                      if(delTPos<maxResid   &&
                         delSlope< maxSlopeDiff && 
                         track1->GetLast()-nremove1>3 && track2->GetLast()-nremove2>3){
                        merge=true;
                        MSG("TrackSR",Msg::kDebug) << " merging tracks " << endl;
                      }         
                    }
                    if(merge){
                      if( nremove1 || nremove2 )  MSG("AlgTrackSRList",Msg::kDebug) << " removing hits " << nremove1 << " " << nremove2 <<  endl;
                      for(Int_t i=0;i<nremove1;i++){
                        if(direction>0){
                          track1->RemoveAt(track1->GetLast());
                        }
                        else{
                          track1->RemoveAt(0);
                        }
                        track1->Compress();
                      }
                      for(Int_t i=0;i<nremove2;i++){
                        if(direction>0){
                          track2->RemoveAt(0);
                        }
                        else{
                          track2->RemoveAt(track2->GetLast());
                        }
                        track2->Compress();
                      }
                    }     
                  }
                  if(merge){
                    for (Int_t icluster=0;icluster<=track2->GetLast();icluster++){
                      TrackClusterSR * tc0=track2->GetCluster(icluster);
                      Double_t slope = track2->GetForwardSlope(icluster); 
                      track1->Add(tc0,slope);
                    }
                    track2->SetIsBad(true);
                  }
                }
              }
            }
            itrack2++;
          }
        }
        itrack1++;
      }
    }
  }
    //remove the bad tracks from the list
  for (Int_t itrk=0; itrk<=tracklist->GetLast(); itrk++) {
    track1 = dynamic_cast<Track2DSR*>(tracklist->At(itrk));
    if(track1->IsBad()){
      MSG("TrackSR" , Msg::kDebug) << "removing bad track " << itrk 
                                   << " with beg plane = "
                                   << track1->GetBegPlane() << endl;
      
      tracklist->RemoveAt(itrk);
      delete track1;
    }
  }
  tracklist->Compress();
}


//-----------------------------------------------------------------------
//this method fills in the gaps that may remain in the 2d tracks
void AlgTrackSRList::FillInGaps(Track2DSRItr &trackItr, 
                                TrackClusterSRItr &clusterItr,
                                Int_t direction)
{
  Track2DSR *thistrack = 0;
  MSG("TrackSR",Msg::kDebug) << "filling in gaps" << endl;

  //array of clusters to add to the track
  TObjArray *addtotrack = new TObjArray(1,0);

  TrackClusterSR *tc = 0;
  TrackClusterSR *tc0 = 0;
  TrackClusterSR *tc1 = 0;
  TrackClusterSR *besttc = 0;

  clusterItr.ResetFirst();
  
  PlaneView::PlaneView_t view = PlaneView::kUnknown;
  if(clusterItr.Ptr() ) view = clusterItr.Ptr()->GetPlaneView();

  Int_t numAdded = 0;

  Double_t dslope = 0.;
  Double_t tcerr0 = 0.;
  Double_t tcerr1 = 0.;
  Double_t tcerr = 0.;
  Double_t bestresid = 99999999.;
  Double_t proj = 0.;
  Double_t dtpos = 0.;
 
  trackItr.ResetFirst();
  clusterItr.Reset();
  

  if(direction == 1)clusterItr.ResetFirst();
  else if(direction == -1) clusterItr.ResetLast(); 
  trackItr.ResetFirst();

  //loop over tracks
  while( (thistrack = trackItr()) ){

    //on a new track so clear the addtotrack array
    addtotrack->Clear();
    numAdded = 1;

    map<TrackClusterSR*,Double_t> thisslope;

    while( numAdded>0 && thistrack->GetPlaneView() == view){
      
      //reset numAdded to 0 for the loop
      numAdded = 0;
      
 
        
      MSG("TrackSR",Msg::kDebug) << "track beg plane = " 
                                        << thistrack->GetBegPlane() 
                                        << " end plane = " 
                                        << thistrack->GetEndPlane() 
                                 << endl;
      
      //loop over the clusters in this track
      for(int ihit=0; ihit<thistrack->GetLast(); ++ihit){
        tc0 = thistrack->GetCluster(ihit);
        tc1 = thistrack->GetCluster(ihit+1);
        
        //make sure clusters are not from consecutive planes in the same view
        //also if one of the clusters is a track end point, look for clusters
        //beyond it
        if(TMath::Abs(tc0->GetPlane()-tc1->GetPlane())>2
           || ihit == 0 || ihit+1 == thistrack->GetLast() ){
          
          MSG("TrackSR",Msg::kDebug) << "tc0 " << tc0->GetPlane() << " " 
                                            << tc0->GetMinStrip() << " " 
                                            << tc0->GetMaxStrip() << " " 
                                            << tc0->GetTPos() << " " 
                                            << tc0->GetTPosError() << endl 
                                            << "tc1 " << tc1->GetPlane() << " " 
                                            << tc1->GetMinStrip() << " " 
                                            << tc1->GetMaxStrip() << " " 
                                            << tc1->GetTPos() << " " 
                                            << tc1->GetTPosError() << endl;
          
          dslope = (tc0->GetTPos()-tc1->GetTPos())/(tc0->GetZPos()-tc1->GetZPos());
          tcerr0 = tc0->GetTPosError();
          tcerr1 = tc1->GetTPosError();
          tcerr = sqrt(tcerr0*tcerr0+tcerr1*tcerr1);
          
          // limit tcerr to 5 cm
          if (tcerr>0.05) tcerr = 0.05;
          
          MSG("TrackSR",Msg::kDebug) << "dslope " << dslope << "  errors = " 
                                     << tcerr0 << " " << tcerr1 << " " << tcerr 
                                     << endl;
          
          //on a new gap, so reset the best residual and cluster variables
          bestresid = 99999.;
          besttc = 0;
          
          //loop over the master list of track clusters
          if(direction == 1)clusterItr.ResetFirst();
          else if(direction == -1) clusterItr.ResetLast();
          
          while( (tc = clusterItr()) ){

            
            //check to see if tc comes from the plane before the first in 
            //the track and isnt too far away from the first hit plane
            
            //or if tc comes from the plane after the last in
            //the track and isnt too far away from the last hit plane
            
            //or if tc comes from a plane between tc0 and tc1
            
            //remember, the track's clusters are in order according to the
            //direction, so the plane of tc0 is before that of tc1 in time


            if((tc->GetPlane()*direction>tc0->GetPlane()*direction
                && tc->GetPlane()*direction<tc1->GetPlane()*direction)
                 || (ihit+1 == thistrack->GetLast() 
                     && tc->GetPlane()*direction>tc1->GetPlane()*direction
                     && TMath::Abs(tc->GetPlane()
                                   -tc1->GetPlane())<=fParmTrk2DNSkip*TMath::Abs(tc0->GetPlane()-tc1->GetPlane()))
                 || (ihit == 0 && 
                     tc->GetPlane()*direction<tc0->GetPlane()*direction
                     && TMath::Abs(tc->GetPlane()
                                   -tc0->GetPlane())<=fParmTrk2DNSkip*TMath::Abs(tc0->GetPlane()-tc1->GetPlane()))){
              
                
                //get the projected location for the plane of tc from the slope
                //between the two bounding clusters
                
                proj = tc0->GetTPos()+dslope*(tc->GetZPos()-tc0->GetZPos());
                
                dtpos = 0.;
                
                //if the transverse position of tc is out of bounds based on the
                //projected location +- the error, then set dtpos to the minimum
                //of the difference in its tpos +- its error - projection +- error
                if(tc->GetTPos()<proj-tcerr || tc->GetTPos()>proj+tcerr)
                  dtpos=min(TMath::Abs(tc->GetTPos()+tc->GetTPosError()
                                       -(proj+tcerr)),
                            TMath::Abs(tc->GetTPos()-tc->GetTPosError()
                                       -(proj-tcerr)));
                
                MSG("TrackSR",Msg::kDebug) << "considering " 
                                                  << tc->GetPlane() 
                                                  << " tpos " << tc->GetTPos() 
                                                  << " proj " << proj 
                                                  << " dtpos " << dtpos 
                                           << "/" <<  fParmTrk2DMaxResid*(1+0.15*TMath::Abs(tc->GetPlane()-tc0->GetPlane()))
                                                  << endl;
                
                //if the dtpos is acceptable and less than the best residual so far
                Float_t dplaneproj=min(TMath::Abs(tc->GetPlane()-tc0->GetPlane()),TMath::Abs(tc->GetPlane()-tc1->GetPlane()));
                if(dtpos<fParmTrk2DMaxResid*(1+0.15*dplaneproj) && dtpos<bestresid){
                  besttc = tc;
                  bestresid = dtpos;
                }
                //            }
            }//end if on a plane between the bounding clusters
          }//end loop over master cluster list
          
          //you found a cluster to add, so add it
          if(besttc){
            addtotrack->AddLast(besttc);
            thisslope[besttc] = ((tc0->GetTPos()-tc1->GetTPos())
                                 /(tc0->GetZPos()-tc1->GetZPos()));
            ++numAdded;
            MSG("TrackSR",Msg::kDebug) << "adding " << besttc->GetPlane() 
                                              << " " << besttc->GetTPos() 
                                              << " numAdded = " << numAdded <<endl;
          }
        }//end if in a gap between clusters
      }//end loop over clusters in the track

      
      //now that you have found all the clusters to add to the track, add them
      for(int i=0; i<=addtotrack->GetLast(); ++i) {
        tc1 = dynamic_cast<TrackClusterSR*>(addtotrack->At(i));
        MSG("TrackSR",Msg::kDebug) << " adding " << tc1->GetPlane() 
                                          << " " << tc1->GetTPos() << endl;
        thistrack->Add(tc1,thisslope[tc1]);

        //added the cluster, now remove it from the addtotrack array
        addtotrack->RemoveAt(i);
      }

      addtotrack->Compress();
    
      }//end loop to add clusters to the track


      // remove isloated clusters at track ends i
    if(thistrack->GetLast()>2){
      Float_t dz_front=TMath::Abs(thistrack->GetCluster(0)->GetZPos()-
                            thistrack->GetCluster(1)->GetZPos());
      while(dz_front>3.0 && thistrack->GetLast()>2){
        thistrack->RemoveAt(0);
        dz_front=TMath::Abs(thistrack->GetCluster(0)->GetZPos()-
                      thistrack->GetCluster(1)->GetZPos());
      }
    }
    if(thistrack->GetLast()>2){
      Float_t dz_back=TMath::Abs(thistrack->GetCluster(thistrack->GetLast())->GetZPos()-
                           thistrack->GetCluster(thistrack->GetLast()-1)->GetZPos());
      while(dz_back>3.0 && thistrack->GetLast()>2){
        thistrack->RemoveAt(thistrack->GetLast());
        dz_back=TMath::Abs(thistrack->GetCluster(thistrack->GetLast())->GetZPos()-
                     thistrack->GetCluster(thistrack->GetLast()-1)->GetZPos());    
      }
    }
    Float_t totcharge = 0.;
    for (int i=0; i<=thistrack->GetLast(); i++) {
      tc1 = thistrack->GetCluster(i);
      MSG("TrackSR",Msg::kDebug) << " " << tc1->GetPlane() 
                                        << " " << tc1->GetTPos() << endl;
      totcharge += tc1->GetCharge();
    }
    //     MSG("TrackSR",Msg::kVerbose) << "total charge = " << totcharge << "\n";
    
    }//end loop over tracks

  delete addtotrack;
  return;
}

 
//---------------------------------------------------------------
//this method forms CandTrackSR's out of Track2DSR's.  it assumes
//that all bad tracks have been removed from the list of tracks
void AlgTrackSRList::FormCandTrackSR(Track2DSRItr &uTrackItr, 
                                     Track2DSRItr &vTrackItr,
                                     TObjArray *candTracks, CandHandle &ch,
                                     CandSliceHandle *slice, AlgHandle ah,
                                     CandContext cxx, Int_t direction)
{
  Track2DSR *tracku = 0;
  Track2DSR *trackv = 0;

  uTrackItr.ResetFirst();
  vTrackItr.ResetFirst();
  Int_t nUtracks=0;
  Int_t nVtracks=0;
  while( (tracku = uTrackItr()) )nUtracks++;  
  while( (trackv = vTrackItr()) )nVtracks++;
  uTrackItr.ResetFirst();
  vTrackItr.ResetFirst();

  TrackClusterSR *tc = 0;

  Double_t nhit0;
  Bool_t match = false;
  Int_t nhitoverlap = 0;
  Int_t uTrackCtr = 0, vTrackCtr = 0;

  //loop over the tracks from the 2 views


  while( (tracku = uTrackItr()) ){

   MSG("TrackSR",Msg::kDebug) << "forming 3D tracks, u track " 
                                     << uTrackCtr << " beg plane = " 
                                     << tracku->GetBegPlane()
                                     << " end plane = " << tracku->GetEndPlane() 
                                     << endl;


   vTrackCtr = 0;

    while( (trackv = vTrackItr()) ){

      match=false;
      Int_t begMisses=0;
      Int_t endMisses=0;
      for(Int_t iu=0; iu<=tracku->GetLast(); ++iu){

        tc = tracku->GetCluster(iu);
        if(direction*tc->GetPlane()<direction*trackv->GetBegPlane())begMisses++;
        if(direction*tc->GetPlane()>direction*trackv->GetEndPlane())endMisses++; 
      }//end loop over u clusters
      
      for(Int_t iv=0; iv<=trackv->GetLast(); ++iv){
        tc = trackv->GetCluster(iv);
        if(direction*tc->GetPlane()<direction*tracku->GetBegPlane())begMisses++;
        if(direction*tc->GetPlane()>direction*tracku->GetEndPlane())endMisses++; 
      }//end loop over v clusters
      
      
      MSG("TrackSR",Msg::kDebug) << " # skip Beg/End " << begMisses << " " << endMisses << endl;

      MSG("TrackSR",Msg::kDebug)  << " u track " << uTrackCtr << "/" << nUtracks  << " v track " << vTrackCtr << "/" << nVtracks  
                                        << " beg plane " << trackv->GetBegPlane() 
                                       << " end plane " << trackv->GetEndPlane() 
                                       << endl
                                       << " dtime "
                                       << TMath::Abs(tracku->GetT0()
                                                     -trackv->GetT0()) 
                                       << "/" 
                                       << fParmDiffViewTimeMatch << " beg Misses" 
                                       << begMisses
                                       << "/" << fParmDiffViewBegPlaneMatch 
                                       << " end Misses  " 
                                       << endMisses 
                                       << "/" << fParmDiffViewEndPlaneMatch <<endl;
      
      // check 2D tracks for position and time agreement. If there
      // is only 1 2D track in each view, skip time agreement check
      if((TMath::Abs(tracku->GetT0()-trackv->GetT0())
          <=fParmDiffViewTimeMatch || 
          (nUtracks==1 && nVtracks==1))  
         &&(begMisses<=fParmDiffViewBegPlaneMatch 
            ||  endMisses <= fParmDiffViewEndPlaneMatch)
         && (TMath::Abs(tracku->GetBegPlane()-trackv->GetBegPlane())<= fParmDiffViewBegPlaneMatch*6 || 
             TMath::Abs(tracku->GetEndPlane()-trackv->GetEndPlane())<= fParmDiffViewEndPlaneMatch*6)
         ){
                 
        // check hit overlap
        nhitoverlap = 0;
        //set nhit0 to the minimum of the number of clusters (ie planes) 
        //from each view times the amount of overlap required
        nhit0 = (1.*TMath::Min(tracku->GetLast()+1,trackv->GetLast()+1)
                 *fParmDiffViewPlaneOverlap);
        //loop over the clusters in the u track. count up the number of 
        //clusters in the u view that are between the first and last planes 
        //in the v view track
        for(Int_t iu=0; iu<=tracku->GetLast() && 1.*nhitoverlap<nhit0; ++iu){
          tc = tracku->GetCluster(iu);
          if(direction*tc->GetPlane()>=direction*trackv->GetBegPlane()
             && direction*tc->GetPlane()<=direction*trackv->GetEndPlane()) 
            ++nhitoverlap;
        }//end loop over u clusters
        
        if(nhitoverlap>=nhit0) match = true;

        if(!match){
          nhitoverlap = 0;
          //no match found, loop over v clusters
          for(Int_t iv=0; iv<=trackv->GetLast() && 1.*nhitoverlap<nhit0; ++iv){
            tc = trackv->GetCluster(iv);
            if(direction*tc->GetPlane()>=direction*tracku->GetBegPlane()
               && direction*tc->GetPlane()<=direction*tracku->GetEndPlane()) 
              ++nhitoverlap;
          }//end loop over v clusters

          if (nhitoverlap>=nhit0) match = true;

        }//end if !match
       

        //does it match now?, if so, make a CandTrackSRHandle
        if(match){
          TObjArray trackin;
          trackin.Add(tracku);
          trackin.Add(trackv);
          trackin.Add(slice);
          cxx.SetDataIn(&trackin);
          CandTrackSRHandle trackhandle = CandTrackSR::MakeCandidate(ah,cxx);
          trackhandle.SetCandSlice(slice);
          // count number of track hits with 3D location outside detector
          UgliGeomHandle ugh(fVldc);
          CandStripHandleItr stripItr(trackhandle.GetDaughterIterator());
          Float_t nstrip=0;
          Float_t noutside=0;
          CandStripHandle *strip=0;
          while( (strip = stripItr()) ){
            PlexStripEndId stripid = strip->GetStripEndId();
            UgliStripHandle striphandle = ugh.GetStripHandle(stripid);
            Double_t halflength = striphandle.GetHalfLength();
            Int_t plane=strip->GetPlane();
            Double_t stripX=0;
            switch (strip->GetPlaneView()) {
            case PlaneView::kX: case PlaneView::kU:
              stripX = - trackhandle.GetV(plane);
              break;
            case PlaneView::kY: case PlaneView::kV:
              stripX = trackhandle.GetU(plane);
              break;
            default:
              continue;
            }
            if(fDetector==Detector::kNear) {
              const double rsqrt2 = 1.0/sqrt(2.0);
              const TVector3 kU(rsqrt2,rsqrt2,0.);
              const TVector3 kV(-rsqrt2,rsqrt2,0.);  
              if (strip->GetPlaneView()==PlaneView::kU) {
                double lambda = kV.Dot(striphandle.GlobalPos(0.));
                stripX += lambda;
              } else {
                double lambda = kU.Dot(striphandle.GlobalPos(0.));
                stripX -= lambda;
              }
            }
            if(TMath::Abs(stripX)>halflength+1.0*Munits::m){
              noutside++;
            }
            nstrip++;
          }
          
          // if less than 3/4 of hit planes are within active detector volume, don't keep this track
           if(noutside/nstrip<0.25){
            // the following should be cast to CandTrackSRHandle, but doing 
            // so gives a compilation error
            const CandTrackSRHandle *track = dynamic_cast<const CandTrackSRHandle*>
              (ch.AddDaughterLink(trackhandle));
            candTracks->AddLast(const_cast<CandTrackSRHandle*>(track));
              }
           else{
             //      delete trackhandle;
             MSG("TrackSR", Msg::kDebug) << " **** removing track with " << noutside/nstrip << " outside detector " << endl;
           }
         
        }
//      else {
//        MSG("TrackSR",Msg::kDebug) << "no match, plane overlap = " 
//                                   << nhitoverlap 
//                                   << "/" << nhit0 << "\n";
//      }
      }//end if tracks could match based on extent
      ++vTrackCtr;
    }//end loop over v tracks
   
    //reset the v track iterator
    vTrackItr.ResetFirst();

    ++uTrackCtr;
  }//end loop over u tracks

  MSG("TrackSR", Msg::kDebug) << candTracks->GetLast()+1 
                              << " tracks formed" << endl;

  return;
}

//----------------------------------------------------------------------
//this method deletes twin tracks - those that may be duplicates
void AlgTrackSRList::DeleteTwinTracks(CandTrackSRHandleItr &candTracks1,
                                      CandTrackSRHandleItr &candTracks2,
                                      TObjArray *candTracks, Int_t uMaxTrack,
                                      Int_t vMaxTrack, Int_t nmax3dtrk,
                                      CandHandle &ch)
{

//   MSG("TrackSR",Msg::kDebug) << "deleting twin tracks" << endl;;
 

  Int_t ncompat=0;
  Bool_t * compatlist = new Bool_t [candTracks1.SizeSelect()];

  //initialize the array
  Int_t ntrk = candTracks1.SizeSelect();

  for(int i=0;i<ntrk;i++){
    compatlist[i]=false;
  }
  
  candTracks1.ResetFirst();
  candTracks2.ResetFirst();

  CandTrackSRHandle *track0 = 0;
  CandTrackSRHandle *besttrack = 0;
  
  Int_t itrk0 = 0;

  Int_t besttrk = 0;
  Int_t bestlen=0;
  Float_t bestchi2=1e9;
  Int_t bestnplanes=0;
  // find the best track, based on len, and chi2
  while( (track0 = candTracks1()) ){
    Int_t nplanes=0;
    Int_t len=0;
    Float_t chi2=0;   
    len =TMath::Abs(track0->GetBegPlane()-track0->GetEndPlane());
    nplanes = track0->GetNPlane();
    chi2 =track0->GetUTrack()->GetChi2()+track0->GetVTrack()->GetChi2();
    if(nplanes>bestnplanes){
      bestnplanes=nplanes;
      bestlen=len;
      bestchi2=chi2;
      besttrack=track0;
      besttrk=itrk0;
    }
    else if (nplanes==bestnplanes){
      if(len>bestlen){
        bestnplanes=nplanes;
        bestlen=len;
        bestchi2=chi2;
        besttrack=track0;
        besttrk=itrk0;
      }
      else if( len==bestlen){
        if(chi2<bestchi2){
          bestnplanes=nplanes;
          bestlen=len;
          bestchi2=chi2;
          besttrack=track0;
          besttrk=itrk0;
        }
      }
    }
    itrk0++;
  }
  if(besttrack){
    MSG("TrackSR", Msg::kDebug) << "best track  " << besttrk << " extent: " << besttrack->GetEndPlane() << " " << besttrack->GetBegPlane() << endl;
  }
  if(!besttrack)return;
  ncompat=1;
  compatlist[besttrk]=true;

  // now find all other tracks compatible with the best
  candTracks1.Reset();
  itrk0=0;
  Int_t itrk1=0;
  const CandTrackSRHandle * track1=0;
  while( (track0 = candTracks1()) ){
    MSG("TrackSR", Msg::kDebug) << "outer loop track " << itrk0 << endl;
    Bool_t trackcompat=true;    
    for(int itrk=0;itrk<ntrk;itrk++){
      itrk1=itrk;
      TObject *tobj1 = candTracks->At(itrk);
      track1 = dynamic_cast<CandTrackSRHandle*>(tobj1);
      if(compatlist[itrk1] && *track0!=*track1){
        MSG("TrackSR", Msg::kDebug) << "inner loop track " << itrk1 << endl;
        Double_t x0 = 0.;
        Double_t x1 = 0.;
        Double_t x01 = 0.;
        
        Double_t x0v = 0.;
        Double_t x1v = 0.;
        Double_t x01v = 0.;
        
        Double_t x0u = 0.;
        Double_t x1u = 0.;
        Double_t x01u = 0.;
        x0u=track0->GetUTrack()->GetLast()+1;
        x1u=track1->GetUTrack()->GetLast()+1;
        TrackClusterSR *clus1;
        TrackClusterSR *clus2;
        for(Int_t i = 0; i<=track0->GetUTrack()->GetLast(); ++i){
          clus1 = track0->GetUTrack()->GetCluster(i);
          for(Int_t j = 0; j<=track1->GetUTrack()->GetLast(); ++j){
            clus2 = track1->GetUTrack()->GetCluster(j);
            if((clus1->GetPlane()==clus2->GetPlane()) && (clus1->GetMinStrip()==clus2->GetMinStrip()) && 
               (clus1->GetMaxStrip()==clus2->GetMaxStrip()) && (clus1->GetNStrip()==clus2->GetNStrip())) {
              x01u++;
            }
          }      
        }        
        
        x0v=track0->GetVTrack()->GetLast()+1;
        x1v=track1->GetVTrack()->GetLast()+1;
        for(Int_t i = 0; i<=track0->GetVTrack()->GetLast(); ++i){
          clus1 = track0->GetVTrack()->GetCluster(i);
          for(Int_t j = 0; j<=track1->GetVTrack()->GetLast(); ++j){
            clus2 = track1->GetVTrack()->GetCluster(j);
            if((clus1->GetPlane()==clus2->GetPlane()) && (clus1->GetMinStrip()==clus2->GetMinStrip()) && 
               (clus1->GetMaxStrip()==clus2->GetMaxStrip()) && (clus1->GetNStrip()==clus2->GetNStrip())) {              
              x01v++;
            }
          }      
        }   
        x0 = 1.*track0->GetNStrip();
        x1 = 1.*track1->GetNStrip();
        x01 = 1.*track0->GetNStrip(track1);                       
        
        MSG("TrackSR", Msg::kDebug) << " V overlap " << x01v/x0v << "/" << fParmTrk3DTwinFrac << "/" << x01v/x1v << " U overlap " << x01u/x0u << "/" << fParmTrk3DTwinFrac << "/" << x01u/x1u << endl;  
        MSG("TrackSR", Msg::kDebug) << " 3d  overlap " << x01/x0 << "/" << fParmTrk3DTwinFrac << "/" << x01/x1  << endl;
        if(x01v/x0v>=fParmTrk3DTwinFrac ||
           x01v/x1v>=fParmTrk3DTwinFrac ||
           x01u/x0u>=fParmTrk3DTwinFrac ||
           x01u/x1u>=fParmTrk3DTwinFrac ||
           x01/x0>=fParmTrk3DTwinFrac  ||
           x01/x1>=fParmTrk3DTwinFrac) {
          MSG("TrackSR", Msg::kDebug) << " tracks have too much overlap - eliminate track  " << itrk0 << endl;
          trackcompat=false;
          Int_t len0 =TMath::Abs(track0->GetBegPlane()-track0->GetEndPlane());
          Int_t nplanes0 = track0->GetNPlane();
          Float_t chi2_0 =track0->GetUTrack()->GetChi2()+track0->GetVTrack()->GetChi2();
          
          Int_t len1 = TMath::Abs(track1->GetBegPlane()-track1->GetEndPlane());
          Int_t nplanes1 = track1->GetNPlane();
          Float_t chi2_1 =track1->GetUTrack()->GetChi2()+track1->GetVTrack()->GetChi2();
          
          if(nplanes0>nplanes1){
            trackcompat=true;
            compatlist[itrk0]=true;
            compatlist[itrk1]=false;
            MSG("TrackSR", Msg::kDebug) << " eliminating track  " << itrk1 << " instead " << endl;
          }
          else if (nplanes0==nplanes1){
            if(len0>len1){
              trackcompat=true;
              compatlist[itrk0]=true;
              compatlist[itrk1]=false;
              MSG("TrackSR", Msg::kDebug) << " eliminating track  " << itrk1 << "instead " << endl;
            }
            else if( len0==len1){
              if(chi2_0<chi2_1){
                trackcompat=true;
                compatlist[itrk0]=true;
                compatlist[itrk1]=false;
                MSG("TrackSR", Msg::kDebug) << " eliminating track " << itrk1 << "instead " << endl;
              }
            }
          }
          if(!compatlist[itrk0])break;
        }
      }
    }
    if(trackcompat && !compatlist[itrk0])ncompat++;
    compatlist[itrk0]=trackcompat;
    itrk0++;
  }
  MSG("TrackSR", Msg::kDebug) << ncompat << " compatible tracks " << endl;
  // sort tracks by strip multiplicity
  Int_t *iassigned = new Int_t[candTracks->GetLast()+1];
  for (Int_t indx=0; indx<=candTracks->GetLast(); ++indx) {
    iassigned[indx] = 0;
  }

  Int_t nmax3d = 0;
  if(fParmTrk3DTrackMax>0) nmax3d = fParmTrk3DTrackMax;
  else{
    nmax3d = max(uMaxTrack,vMaxTrack);
    nmax3d = max(nmax3d,nmax3dtrk);
  }

  // allow for at least one track
  if (nmax3d<1) nmax3d = 1;

  MSG("TrackSR", Msg::kDebug) << "nmax3d = " << nmax3d <<  " p3DTrackMax  " << fParmTrk3DTrackMax << endl;

  //loop over all the tracks to find the n tracks with the most
  //hit strips where n = minimum(nmax3d, candTracks->GetLast()+1
  Int_t nLoop = min(nmax3d,candTracks->GetLast()+1);
  for(Int_t i = 0; i < nLoop; ++i){
    int imax=-1;
    Int_t nmax = 0;
    
    //loop over the tracks that werent removed to find which one has 
    //the largest number of strips - that will be the one we eventually keep
    for (itrk0=0; itrk0<=candTracks->GetLast(); ++itrk0){
      TObject *tobj0 = candTracks->At(itrk0);
      const CandTrackSRHandle *track0 =
        dynamic_cast<CandTrackSRHandle*>(tobj0);
      if (compatlist[itrk0] && !iassigned[itrk0]){
        if (track0->GetNStrip()>nmax) {
          nmax = track0->GetNStrip();
          imax = itrk0;
        }
      }
    }//end loop over tracks
    if(imax>=0) iassigned[imax] = 1;
  }//end loop to get n tracks with most hit strips

  //get rid of all tracks but the one with the most strips hit in it
  for(itrk0=0; itrk0<=candTracks->GetLast(); ++itrk0){
    if(!iassigned[itrk0]){
      // remove from candidate daughter list
      MSG("TrackSR", Msg::kDebug) << "removing track at position " << itrk0+1 
                                  << " of "
                                 << candTracks->GetLast()+1 << endl;
      ch.RemoveDaughter(dynamic_cast<CandTrackSRHandle*>(candTracks->At(itrk0)));
    }
  }
  
  delete [] iassigned;
  delete [] compatlist;
  return;
}

//-----------------------------------------------------------------------
//this method cleans out the intermediate objects used in forming the 
//tracks - the Track2DSR and TrackClusterSR objects
void AlgTrackSRList::CleanUp(TObjArray *trackList, TObjArray *clusterList)
{
  Track2DSR *track2d = 0;

  for(Int_t i=0; i<=trackList->GetLast(); ++i){
    track2d = dynamic_cast<Track2DSR*>(trackList->At(i));
    trackList->RemoveAt(i);
    delete track2d;
  }
  trackList->Compress();
  
  // delete track clusters
  TrackClusterSR *tc = 0;
  for(Int_t i=0; i<=clusterList->GetLast(); ++i){
   
    tc = dynamic_cast<TrackClusterSR*>(clusterList->At(i));
    fMapIsWide[tc->GetPlane()] = 0;
   
    MSG("TrackSR", Msg::kDebug) << "cleaning up, fMapIsWide for plane " 
                               << tc->GetPlane() << " is " 
                               << fMapIsWide[tc->GetPlane()]
                               << endl;
    delete tc;
  }
  clusterList->Compress();
  
  return;
}

//----------------------------------------------------------------------
//this method loops over an array of strips and looks for the planes
//representing the highest and lowest numbers, as well as the total
//number of different planes represented
void AlgTrackSRList::Find2DTrackEndPoints(CandStripHandleItr &stripItr, 
                                          Int_t &vtxPlane, Int_t &endPlane,
                                          Int_t &numPlanes)
{
  //prevPlane keeps track of which plane you were last on
  Int_t prevPlane = 0;
  CandStripHandle *strip = 0;
  Int_t currentPlane = 0;

  while( stripItr.IsValid() ){
    
    strip = stripItr.Ptr();
    currentPlane = strip->GetPlane();

    MSG("TrackSR", Msg::kDebug) << "plane " << currentPlane << " tpos " 
                               << strip->GetTPos() << " zpos " << strip->GetZPos()
                               << " time " << strip->GetTime() << endl;

    if( prevPlane != currentPlane && 
        strip->GetCharge()>fParmMinStripPulseHeight){
      ++numPlanes;
      if( currentPlane < vtxPlane) vtxPlane = currentPlane;
      if( currentPlane > endPlane) endPlane = currentPlane;
      prevPlane = currentPlane;
    }//end if this is a new plane
    
    //advance the iterator
    stripItr.Next();
  }//end loop over u strips
  
  stripItr.ResetFirst();

  return;
}

//--------------------------------------------------------------------
//
// create TrackClusterSR objects
//
void AlgTrackSRList::MakeTrackClusters(TObjArray *trackClusterList, 
                                       CandStripHandleItr &stripItr,
                                       CandTrackSRListHandle &cth,
                                       Int_t expectedStrips,
                                       Int_t direction)
{
  //array of flags to mark showerlike planes - those with at least one
  //wide cluster
  Int_t showerLike[500];
  for(Int_t i=0; i<500; ++i){
    showerLike[i] = 0;
  }

  MSG("AlgTrackSRList", Msg::kDebug) << "expectedStrips " << expectedStrips
                                     << " direction " << direction << endl;

  Int_t oldplane=0;
  Int_t oldstrip=0;
  Double_t oldtime=0.;
  TrackClusterSR *tc=0;
  stripItr.SetDirection(direction);
  stripItr.Reset();

  //loop over the strips in the Iterator
  //declare a pointer and zero it so you dont have to keep declaring 
  //the strip handle
  CandStripHandle *strip = 0;

  while((strip = stripItr())){

    // don't include spectrometer hits
    if(fDetector!=Detector::kNear || strip->GetPlane()<121){
      MSG("TrackSR", Msg::kDebug) << "looking at plane " << strip->GetPlane()
                                  << " strip " << strip->GetStrip() << endl; 
      
      //check to see if you have a current cluster 
      if(tc){
        //now see if this strip comes from the same plane and 
        //is part of the current cluster (ie not separated by too many strips
        if(strip->GetPlane()==oldplane &&
           direction*strip->GetStrip()<=direction*oldstrip+1+fParmTrkClsNSkip){
          tc->AddStrip(strip);
          oldstrip = strip->GetStrip();
          oldtime = strip->GetBegTime();
          
          MSG("TrackSR", Msg::kDebug) << "add strip " << oldstrip 
                                      << " with time "  << oldtime*1.e9 
                                      << " and charge " << strip->GetCharge() 
                                      << " to cluster on plane " 
                                      << oldplane << " " << fParmTrkClsNSkip 
                                      << endl;
          
        }
        else{
          //you have a new plane or cluster - see if the current
          //cluster passes the pusle height cut before adding it to the array
          MSG("TrackSR", Msg::kDebug) << "checking tc charge " << tc->GetCharge() << "/" << fParmMinClusterPulseHeight << endl;
        

          if(tc->GetCharge()>=fParmMinClusterPulseHeight){
            
            //check to see if this is a wide cluster - ie it has more strips
            //in it than expected based on the hough transform - the number of 
            //expected strips is passed into the method.  if it is wide, set the 
            //cluster's wide flag. also set the flag in the showerLike array
            //for this plane
            if(tc->GetNStrip()>fParmMaxNExtraStrip+expectedStrips){
              tc->SetIsWide(true);
              showerLike[tc->GetPlane()] = 1;
              
              MSG("TrackSR", Msg::kDebug) << "setting cluster on plane " 
                                          << tc->GetPlane() << " wide with "
                                          << tc->GetNStrip() << " strips " 
                                          << "expected " << expectedStrips << endl;
              
            }
            
            MSG("TrackSR", Msg::kDebug) << "adding new cluster on plane " 
                                        << tc->GetPlane() << " to list " 
                                        << tc->GetTPos() << " " 
                                        << tc->GetZPos() << " with charge " << tc->GetCharge()
                                        << endl;
            
            trackClusterList->AddLast(tc);
            cth.AddTrackCluster(tc);
          }//end if enough charge in this cluster to add it to the list
          else{
            delete tc;
          }
          //now make a new cluster
          TrackClusterSR *trackcluster = new TrackClusterSR(strip, 
                                                            fParmMisalignmentError);
          oldplane = strip->GetPlane();
          oldstrip = strip->GetStrip();
          oldtime = strip->GetBegTime();
          tc = trackcluster;
          
          MSG("TrackSR", Msg::kDebug) << "make cluster with strip " << oldstrip 
                                      << " with time " << oldtime*1e9 
                                      << " with charge " << strip->GetCharge() 
                                      << " on plane " << oldplane << endl;
          
        }
      }//end if you have a current cluster
      else {//make a new cluster
        TrackClusterSR *trackcluster = new TrackClusterSR(strip, 
                                                          fParmMisalignmentError);
        oldplane = strip->GetPlane();
        oldstrip = strip->GetStrip();
        oldtime = strip->GetBegTime();
        tc = trackcluster;
        
        MSG("TrackSR", Msg::kDebug) << "make cluster with strip " << oldstrip 
                                    << " with time " << oldtime*1e9  
                                    << " with charge " << strip->GetCharge() 
                                    << " on plane " << oldplane << endl;
        
      }
    }
  }//end loop over strips


  //done with loop, but last cluster needs to be added
  if(tc){
    if(tc->GetCharge()>=fParmMinClusterPulseHeight){
      
      //check to see if this is a wide cluster - ie it has more strips
      //in it than expected based on the hough transform - the number of 
      //expected strips is passed into the method.  if it is wide, set the 
      //cluster's wide flag. also set the flag in the showerLike array
      //for this plane
      if(tc->GetNStrip()>fParmMaxNExtraStrip+expectedStrips){
        tc->SetIsWide(true);
        showerLike[tc->GetPlane()] = 1;
      }
      
      MSG("TrackSR", Msg::kDebug) << "adding new cluster on plane " 
                                  << tc->GetPlane() << " to list " 
                                  << tc->GetTPos() << " " << tc->GetZPos() << endl;
      
      trackClusterList->AddLast(tc);
      cth.AddTrackCluster(tc);
    }
    else{
      delete tc;
    }
  }//end if still a cluster to add to the list

  //loop over the clusters in the array to mark ones in showerlike planes
  for(Int_t j=0; j<trackClusterList->GetLast(); ++j){
    tc = dynamic_cast<TrackClusterSR *>(trackClusterList->At(j)); 
    if(showerLike[tc->GetPlane()] == 1) 
      tc->SetInShowerLikePlane(true);
  }//end loop over clusters 
  
  //debugging output for the cluster list made
   for (Int_t i=0; i<=trackClusterList->GetLast(); ++i) {
     tc = dynamic_cast<TrackClusterSR*>(trackClusterList->At(i));
     MSG("TrackSR",Msg::kDebug) << "trackcluster " << i << " plane " 
                                << tc->GetPlane() << " tpos " << tc->GetTPos() 
                                << " zpos " << tc->GetZPos() << endl
                                << " nstrip " << tc->GetNStrip() << " charge " 
                                << tc->GetCharge() << " wide " 
                                << (Int_t)tc->IsWide()
                                << " showerlike " 
                                << (Int_t)tc->InShowerLikePlane() 
                                << endl;
   }
  
  
  stripItr.ResetFirst();
  return;
}

//______________________________________________________________________
void AlgTrackSRList::Trace(const char * /* c */) const
{
}

---