////////////////////////////////////////////////////////////////////////
// $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
		&