////////////////////////////////////////////////////////////////////////
// $Id: AlgTrackSR.cxx,v 1.71 2007/02/04 06:03:24 rhatcher Exp $
//
// AlgTrackSR.cxx
//
// Begin_Html<img src="../../pedestrians.gif" align=center>
// <a href="../source_warning.html">Warning for beginners</a>.<br> 
//
// AlgTrackSR is a concrete TrackSR Algorithm class.
//
// Author:  R. Lee 2001.02.26
// Revised: B. Rebel 2004.03.31
//
// 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 "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "Algorithm/AlgConfig.h"
#include "Candidate/CandContext.h"
//#include "./AlgTrackSR.h"
#include "CandTrackSR/AlgTrackSR.h"
#include "CandTrackSR/CandTrackSRHandle.h"
#include "CandTrackSR/TrackClusterSR.h"
#include "CandTrackSR/Track2DSR.h"
#include "Conventions/CalTimeType.h"
#include "Conventions/Mphysical.h"
#include "Conventions/Munits.h"
#include "Conventions/PlaneView.h"
#include "DataUtil/GetMomFromRange.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Navigation/XxxItr.h"
#include "Plex/PlexStripEndId.h"
#include "Plex/PlexPlaneId.h"
#include "RecoBase/ArrivalTime.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/LinearFit.h"
#include "RecoBase/PropagationVelocity.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliScintPlnHandle.h"
#include "Validity/VldContext.h"
#include "TMath.h"

NavKey KeyOnStripTPos( const CandStripHandle *cth )
{
  return cth->GetTPos();
}

ClassImp(AlgTrackSR)

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

//______________________________________________________________________
AlgTrackSR::AlgTrackSR()
{
}

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

//______________________________________________________________________
void AlgTrackSR::RunAlg(AlgConfig &ac, CandHandle &ch, CandContext &cx)
{


  MSG("Alg", Msg::kDebug) << "Starting AlgTrackSR::RunAlg()" << endl;

  // get parameters
  fParmNExtraStrip = ac.GetInt("NExtraStrip");
  fParmMinPlanePE = ac.GetDouble("MinPlanePE");
  fParmIsCosmic = ac.GetInt("IsCosmic");
  fParmMinStripPE = ac.GetDouble("MinStripPE");
  fParmMaxClusterSizeTimeFit = ac.GetDouble("MaxClusterSizeTimeFit");
  fParmMaxTimeFitRes = ac.GetDouble("MaxTimeFitRes");
  fParmMinTimeFitSize = ac.GetInt("MinTimeFitSize");
  fParmMaxTimeFitOutlier = ac.GetDouble("MaxTimeFitOutlier");

  assert(ch.InheritsFrom("CandTrackSRHandle"));
  CandTrackSRHandle &cth = dynamic_cast<CandTrackSRHandle &>(ch);

  assert(cx.GetDataIn());

  assert(cx.GetDataIn()->InheritsFrom("TObjArray"));

 
  const TObjArray *tary =
    dynamic_cast<const TObjArray*>(cx.GetDataIn());

  Track2DSR *utrack=0;
  Track2DSR *vtrack=0;
  CandStripHandle *strip = 0;

  Int_t ntrackdigit=0;

  const CandSliceHandle *slice = 0;

  for (Int_t i=0; i<=tary->GetLast(); i++) {
    TObject *tobj = tary->At(i);
    if (tobj->InheritsFrom("CandSliceHandle")) {
      slice = dynamic_cast<CandSliceHandle*>(tobj);
      continue;
    }
    Track2DSR *track = dynamic_cast<Track2DSR*>(tobj);
    assert(track);
    if(!utrack && track->GetPlaneView()==PlaneView::kU) utrack = track;
    if(!vtrack && track->GetPlaneView()==PlaneView::kV) vtrack = track;

  }//end loop to find 2D tracks

  FindStripsInTrack(utrack, vtrack, cth);

  CandStripHandleItr trackStripItr(cth.GetDaughterIterator());
  while( (strip = trackStripItr()) ){
    MSG("AlgTrackSR", Msg::kDebug) << strip->GetPlane() << " " << strip->GetTPos() << endl;
  }
  MSG("AlgTrackSR", Msg::kDebug) << endl;

  assert(slice);
  CandStripHandleItr sliceItr(slice->GetDaughterIterator());
  Int_t maxplane=0;
  while( (strip = sliceItr()) ){
    if(strip->GetPlane()>maxplane) maxplane = strip->GetPlane();
  }
  assert(maxplane);

  //there can be at most 485 planes in either detector - round up
  //to 500
  Float_t planepe[500];
  Float_t zpe[500];
  for (int i=0; i<=maxplane; i++) {
    planepe[i] = 0.;
    zpe[i] = 0.;
  }

  sliceItr.Reset();

  while( (strip = sliceItr()) ){
    if(strip->GetPlane()>=0){
      planepe[strip->GetPlane()] += strip->GetCharge(CalDigitType::kPE);
      zpe[strip->GetPlane()] = strip->GetZPos();
    }
  }

  CandStripHandleItr stripItr(cth.GetDaughterIterator());
  CandStripHandle *firststrip = stripItr.Ptr();

  assert(firststrip);

  const VldContext *vldcontext = firststrip->GetVldContext();
  assert(vldcontext);

  //get the detector type from the VldContext
  fDetector = vldcontext->GetDetector();

  UgliGeomHandle ugh(*vldcontext);

  // calculate transverse positions

  MSG("AlgTrackSR", Msg::kDebug) << "in detector type " 
				 << Detector::AsString(fDetector) 
				 << endl
				 << " u track direction = " 
				 << utrack->GetDirection()
				 << " v track direction = " 
				 << vtrack->GetDirection()
				 << endl;

  if(utrack->GetDirection()>0.) cth.SetTimeSlope(1.);
  else cth.SetTimeSlope(-1.);

  CandRecoHandle *candreco = dynamic_cast<CandRecoHandle *>(&cth);
  assert(candreco);
  
  Int_t begplane = candreco->GetBegPlane();
  Int_t endplane = candreco->GetEndPlane();
  
  cth.SetVtxPlane(begplane);
  cth.SetVtxZ(cth.GetZ(begplane));
  cth.SetVtxU(cth.GetU(begplane));
  cth.SetVtxV(cth.GetV(begplane));
  cth.SetDirCosU(cth.GetDirCosU(begplane));
  cth.SetDirCosV(cth.GetDirCosV(begplane));
  cth.SetDirCosZ(cth.GetDirCosZ(begplane));
  
  cth.SetEndPlane(endplane);
  cth.SetEndZ(cth.GetZ(endplane));
  cth.SetEndU(cth.GetU(endplane));
  cth.SetEndV(cth.GetV(endplane));
  cth.SetEndDirCosU(cth.GetDirCosU(endplane));
  cth.SetEndDirCosV(cth.GetDirCosV(endplane));
  cth.SetEndDirCosZ(cth.GetDirCosZ(endplane));

  SetUVZT(&cth);

  Int_t nshower = 0;
  Int_t ntrackstrip = 0;
  Double_t aveTime=0;
  while( (strip = stripItr()) ){
    nshower = cth.IsInShower(strip);
    aveTime += strip->GetTime(); 
    if(nshower<=1){
      ++ntrackstrip;
      ntrackdigit += strip->GetNDaughters();
    }
  }//end loop to see how many track strips there are

  if(cth.GetNDaughters()>0) aveTime /=cth.GetNDaughters();
  cth.SetVtxT(aveTime);
  cth.SetEndT(aveTime);

  stripItr.Reset();
  
  cth.SetNTrackStrip(ntrackstrip);
  
  //timewalk correction constants
  //  Double_t c1 = ac.GetDouble("TimeWalk1"); //-14.45;
  //  Double_t c2 = ac.GetDouble("TimeWalk2"); //7.498
  //  Double_t c3 = ac.GetDouble("TimeWalk3"); //-1.566;
  Double_t maxPathLength = 0.;
  Double_t timefitchi2 = 0.;
  Double_t fitparm[] = {0.,0.};
 
  
  Int_t ntimedigit = FindTimingDirection(cth, stripItr, utrack, vtrack, 
					 fitparm, maxPathLength, timefitchi2);
 
  if(ntimedigit<=fParmMinTimeFitSize){
    fitparm[0]=aveTime;
    fitparm[1]=0;
    timefitchi2=0;
  }

  stripItr.Reset();

  if(!fParmIsCosmic 
     || fitparm[1]*utrack->GetDirection()>0.)
    cth.SetTimeSlope(TMath::Abs(fitparm[1]));
  else
    cth.SetTimeSlope(-1.*TMath::Abs(fitparm[1]));
  
  if(fitparm[1]<0. || !fParmIsCosmic){ 
    // in beam mode, reconstruct tracks backward, but now want track going forward
    cth.SetTimeOffset(fitparm[0]+maxPathLength*fitparm[1]);
    
    // direction changed from initial guess
    MSG("AlgTrackSR",Msg::kDebug) << "DIRECTION CHANGED from " 
				  << utrack->GetDirection() << endl;
    Int_t vtxplane = cth.GetVtxPlane();
    Int_t endplane = cth.GetEndPlane();
    cth.SetVtxPlane(endplane);
    cth.SetEndPlane(vtxplane);
    
    MSG("AlgTrackSR",Msg::kDebug) << "new vertex and end planes = " 
				  << cth.GetVtxPlane() << " " 
				  << cth.GetEndPlane() << endl;
    SetUVZT(&cth);
  } 
  else
    cth.SetTimeOffset(fitparm[0]);
  
  cth.SetNTrackDigit(ntrackdigit);
  cth.SetNTimeFitDigit(ntimedigit);
  cth.SetTimeFitChi2(timefitchi2);
  
  MSG("AlgTrackSR",Msg::kVerbose) << "timeslope = " << cth.GetTimeSlope() << endl;
  MSG("AlgTrackSR",Msg::kVerbose) << "ntrackdigit = " << cth.GetNTrackDigit() 
				  << endl;
  MSG("AlgTrackSR",Msg::kVerbose) << "ntimefitdigit = " << cth.GetNTimeFitDigit() 
				  << endl;
  MSG("TrackSR",Msg::kVerbose) << "timefitchi2 = " << cth.GetTimeFitChi2() << endl;

  begplane = cth.GetBegPlane();
  endplane = cth.GetEndPlane();

  Track2DSR *utrack0 = utrack->Dup();
  Track2DSR *vtrack0 = vtrack->Dup();

  cth.SetUTrack(utrack0);
  cth.SetVTrack(vtrack0);

  cth.SetVtxPlane(begplane);
  cth.SetVtxZ(cth.GetZ(begplane));
  cth.SetVtxU(cth.GetU(begplane));
  cth.SetVtxV(cth.GetV(begplane));
  cth.SetVtxT(cth.GetTimeOffset());
  cth.SetDirCosU(cth.GetDirCosU(begplane));
  cth.SetDirCosV(cth.GetDirCosV(begplane));
  cth.SetDirCosZ(cth.GetDirCosZ(begplane));

  cth.SetEndPlane(endplane);
  cth.SetEndZ(cth.GetZ(endplane));
  cth.SetEndU(cth.GetU(endplane));
  cth.SetEndV(cth.GetV(endplane));
  cth.SetEndT(cth.GetTimeOffset()+cth.GetTimeSlope()*cth.GetdS());
  cth.SetEndDirCosU(cth.GetDirCosU(endplane));
  cth.SetEndDirCosV(cth.GetDirCosV(endplane));
  cth.SetEndDirCosZ(cth.GetDirCosZ(endplane));

  Int_t plane0 = -999;
  Int_t plane1 = -999;
  Int_t idir = (cth.GetDirCosZ()>0. ? 1 : -1);
  Double_t minTPos = 0.;
  Double_t maxTPos = 0.;
  Double_t minZ    = 0.;
  Double_t maxZ    = 0.;
  ugh.GetTransverseExtent(PlaneView::kU, minTPos, maxTPos);
  ugh.GetZExtent(minZ, maxZ);

  for (int i=0; i<maxplane; i++) {
    if (planepe[i]>=fParmMinPlanePE && planepe[i+1]>=fParmMinPlanePE) {
      if (idir>0) {
        if (plane0<0) {
          plane0 = i;
        }
        if (plane1<0 || i+1>plane1) {
          plane1 = i+1;
        }
      } else {
        if (plane1<0) {
          plane1 = i;
        }
        if (plane0<0 || i+1>plane0) {
          plane0 = i+1;
        }
      }
    }
  }
  MSG("TrackSR",Msg::kDebug) << "Vertex and End Planes >= " << fParmMinPlanePE 
			     << " pe: " << plane0 << " " << plane1 << endl;

  if(fParmIsCosmic){
    Int_t thisplane = cth.GetVtxPlane()-idir;
    Double_t thisu = cth.GetVtxU();
    Double_t thisv = cth.GetVtxV();
    Double_t thisz = cth.GetVtxZ();
    Double_t thisx = 0.707*(thisu-thisv);
    Double_t thisy = 0.707*(thisu+thisv);
    Double_t thisdu = cth.GetDirCosU()/cth.GetDirCosZ();
    Double_t thisdv = cth.GetDirCosV()/cth.GetDirCosZ();
    Double_t thisdx = 0.707*(thisdu-thisdv);
    Double_t thisdy = 0.707*(thisdu+thisdv);
    Int_t nplaneskip = 0;
    while(thisplane>=0 && thisplane<=500 && thisplane*idir>=plane0*idir 
	  && thisx<=maxTPos && thisx>=minTPos 
	  && thisy<=maxTPos && thisy>=minTPos
	  && thisu<=maxTPos && thisu>=minTPos
	  && thisv<=maxTPos && thisv>=minTPos 
	  && nplaneskip<=3
	  ){

      PlexPlaneId plexplane(fDetector,thisplane);
      
      if( plexplane.IsValid() ){
        UgliScintPlnHandle planeid = ugh.GetScintPlnHandle(plexplane);
        if( planeid.IsValid() ){
          Double_t dz = (Double_t)planeid.GetZ0()-thisz;
          if (thisx<=maxTPos && thisx>=minTPos 
	      && thisy<=maxTPos && thisy>=minTPos
	      && thisu<=maxTPos && thisu>=minTPos
	      && thisv<=maxTPos && thisv>=minTPos){ 
	    thisz = (Double_t)planeid.GetZ0();
	    thisx += dz*thisdx;
            thisy += dz*thisdy;
            thisu += dz*thisdu;
            thisv += dz*thisdv;
            cth.SetVtxU(thisu);
            cth.SetVtxV(thisv);
            cth.SetVtxZ(thisz);
            cth.SetVtxPlane(thisplane);
            cth.SetU(thisplane,thisu);
            cth.SetV(thisplane,thisv);
            MSG("TrackSR",Msg::kDebug) << "SETTING VERTEX, (u,v,z) = (" << thisu 
				       << "," << thisv << "," << thisz << ")" << endl;
          }
          if (planepe[thisplane]>0.) {
            nplaneskip = 0;
          } else {
            nplaneskip++;
          }
        }
      }
      thisplane -= idir;
    }
    thisplane = cth.GetEndPlane()+idir;
    thisu = cth.GetEndU();
    thisv = cth.GetEndV();
    thisz = cth.GetEndZ();
    thisx = 0.707*(thisu-thisv);
    thisy = 0.707*(thisu+thisv);
    thisdu = cth.GetEndDirCosU()/cth.GetEndDirCosZ();
    thisdv = cth.GetEndDirCosV()/cth.GetEndDirCosZ();
    thisdx = 0.707*(thisdu-thisdv);
    thisdy = 0.707*(thisdu+thisdv);
    nplaneskip = 0;
    while(thisplane>=0 && thisplane<=500 && thisplane*idir<=plane1*idir 
	  && thisx<=maxTPos && thisx>=minTPos 
	  && thisy<=maxTPos && thisy>=minTPos
	  && thisu<=maxTPos && thisu>=minTPos
	  && thisv<=maxTPos && thisv>=minTPos
	  && nplaneskip<=3
	  ){

      PlexPlaneId plexplane(fDetector,thisplane);
      if( plexplane.IsValid() ){
        UgliScintPlnHandle planeid = ugh.GetScintPlnHandle(plexplane);
        if (planeid.IsValid()) {
          Double_t dz = (Double_t)planeid.GetZ0()-thisz;
          if (thisx<=maxTPos && thisx>=minTPos 
	      && thisy<=maxTPos && thisy>=minTPos
	      && thisu<=maxTPos && thisu>=minTPos
	      && thisv<=maxTPos && thisv>=minTPos){
            thisz = (Double_t)planeid.GetZ0();
            thisx += dz*thisdx;
            thisy += dz*thisdy;
            thisu += dz*thisdu;
            thisv += dz*thisdv;
            cth.SetEndU(thisu);
            cth.SetEndV(thisv);
            cth.SetEndZ(thisz);
            cth.SetEndPlane(thisplane);
            cth.SetU(thisplane,thisu);
            cth.SetV(thisplane,thisv);
            MSG("TrackSR",Msg::kDebug) << "SETTING END, (u,v,z) = (" << thisu 
				       << "," << thisv << "," << thisz << ")" << endl;
          }
          if (planepe[thisplane]>0.) {
            nplaneskip = 0;
          } else {
            nplaneskip++;
          }
        }
      }
      thisplane += idir;
    }
  }


  // traces calculated assuming far detector

  Double_t thistrace[2][4] = {{0.,0.,0.,0.},{0.,0.,0.,0.}}; // u,v,x,y
  Double_t thistracez[2][4] = {{0.,0.,0.,0.},{0.,0.,0.,0.}};

  Double_t thisx[2][4] = {{cth.GetVtxU(),cth.GetVtxV(),0.,0.},
                           {cth.GetEndU(),cth.GetEndV(),0.,0.}};
  Double_t thisz[2] = {cth.GetVtxZ(),cth.GetEndZ()};

  Double_t thisdcos[2][4] = {{-cth.GetDirCosU(),-cth.GetDirCosV(),
                              0.,0.},
                             {cth.GetEndDirCosU(),cth.GetEndDirCosV(),
                              0.,0.}};

  Double_t thisdcosz[2] = {-cth.GetDirCosZ(),cth.GetEndDirCosZ()};

  Double_t mintrace[2] = {999.,999.};
  Double_t mintracez[2] = {999.,999.};

  for (Int_t i=0; i<2; i++) {
    for (Int_t j=0; j<4; j++) {
      if (j==2) {
        thisx[i][j] = .70710678*(thisx[i][0]-thisx[i][1]);
        thisdcos[i][j] = .70710678*(thisdcos[i][0]-thisdcos[i][1]);
      }     
      if (j==3) {
        thisx[i][j] = .70710678*(thisx[i][0]+thisx[i][1]);
        thisdcos[i][j] = .70710678*(thisdcos[i][0]+thisdcos[i][1]);
      }
      if (thisx[i][j]<=maxTPos && thisx[i][j]>=minTPos
	  && thisz[i]>=minZ && thisz[i]<=maxZ) {
        if (thisdcos[i][j]>0.) {
          thistrace[i][j] = (maxTPos-thisx[i][j])/TMath::Abs(thisdcos[i][j]);
        } 
	else if( thisdcos[i][j] < 0.) {
          thistrace[i][j] = (maxTPos+thisx[i][j])/TMath::Abs(thisdcos[i][j]);
	}
        thistracez[i][j] = thistrace[i][j]*TMath::Abs(thisdcosz[i]);
// check in z direction
        Double_t delz = 0.;
        if (thisdcosz[i]>0) {
          delz = maxZ-thisz[i];
        } else {
          delz = thisz[i]-minZ;
        }
        if (delz<thistracez[i][j] && thisdcosz[i] != 0.) {
          thistrace[i][j] = delz/TMath::Abs(thisdcosz[i]);
          thistracez[i][j] = delz;
        }
      }
      if (thistrace[i][j]<mintrace[i]) {
        mintrace[i] = thistrace[i][j];
        mintracez[i] = thistracez[i][j];
      }
    }
  }

  cth.SetVtxTrace(mintrace[0]);
  cth.SetVtxTraceZ(mintracez[0]);
  cth.SetEndTrace(mintrace[1]);
  cth.SetEndTraceZ(mintracez[1]);

  SetdS(&cth);
  Double_t range = cth.GetRange();
  if (range>0.) {
    //    cth.SetMomentum(0.105658*exp(1.0363*log(range)-4.383));
  // taken from PDG R/M vs p/M plot for Fe
  // until June 2006 used to be  cth.SetMomentum(0.048 + 1.660e-3*range + 3.057e-8*range*range);
     cth.SetMomentum(GetMomFromRange(range));   
  } else {
    cth.SetMomentum(0.); // range not valid, set momentum to 0 GeV/c
  }
  Calibrate(&cth);

}

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

//----------------------------------------------------------------------
//figure out which strips belong in the track

void AlgTrackSR::FindStripsInTrack(Track2DSR *uTrack, Track2DSR *vTrack,
				   CandTrackSRHandle &cth)
{
  MSG("AlgTrackSR", Msg::kDebug) << "in find strips in track" << endl;
  
  Track2DSR *track = 0;
  Int_t trkCtr = 0;

  TrackClusterSR *tc = 0;
  CandStripHandle *strip = 0;
  Int_t nstripexp = 0;
  PlaneView::PlaneView_t stripView = PlaneView::kUnknown;
  
  Int_t itcbef=0,itcaft=0;
  TrackClusterSR *tcbef = 0;
  TrackClusterSR *tcaft = 0;
  Double_t localslope = 0.;
  Double_t localtpos = 0.;

  Int_t numStrips = 0;
  
  Float_t stripTPos[192];
  Float_t stripCharge[192];
  Int_t sortedStrips[192];
  
  Int_t indxbest = 0;
  Double_t distbest = 0;
  Double_t dist = 0.;
  Double_t distsum = 0.;
 
  Int_t nshower=0;
  
  while( trkCtr<2 ){
    
    track = uTrack;
    if(trkCtr==1) track = vTrack;

    ++trkCtr;

    //loop over the clusters in the track
    for(Int_t itc=0; itc<track->GetLast()+1; ++itc){
      tc = track->GetCluster(itc);
      numStrips = tc->GetStripList()->GetLast()+1;
      cth.SetTrackPointError(tc->GetPlane(),tc->GetTPosError());
      // assume 4:1 aspect ratio for scintillator strips
      nstripexp = TMath::Max(1,(Int_t)(TMath::Abs(track->GetSlope(itc))/4.)+1);
      nstripexp += fParmNExtraStrip;
      
      //if this cluster has less than the maximum number of strips allowed
      //take all the strips
      if(numStrips<=nstripexp){
	
	MSG("AlgTrackSR", Msg::kDebug) << "cluster on plane " << tc->GetPlane()
				       << " has " 
				       << numStrips
				       << "/" << nstripexp << " strips" << endl;

	//loop over the strips in the cluster
        for(Int_t istrip=0; istrip<numStrips; ++istrip){
          
	  strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(istrip));
	  stripView = strip->GetPlaneView();
	  
          if(stripView == PlaneView::kU || stripView == PlaneView::kV){
	    cth.AddDaughterLink(*strip);
	    cth.SetInShower(strip,0);
	  }
        }
      }
      else{ 
	
	//more than the expected number of strips,
	//only add strips which match fit
	
	MSG("AlgTrackSR", Msg::kDebug) << "cluster on plane " << tc->GetPlane()
					 << " has " 
					 << tc->GetStripList()->GetLast()+1
					 << "/" << nstripexp << " strips" << endl;

	itcbef=0;
	itcaft=0;
	
	//get the bounding clusters.  if on an endpoint cluster get the two
	//closest to it
	if(itc == 0){
	  itcbef = FindNeighborIndex(track, itc, 1);
	  itcaft = FindNeighborIndex(track, itcbef, 1);
	}
	else if(itc == track->GetLast() ){
	  itcbef = FindNeighborIndex(track, itc, -1);
	  itcaft = FindNeighborIndex(track, itcbef, -1);
	}
	else{
	  itcbef = FindNeighborIndex(track, itc,-1);
	  itcaft = FindNeighborIndex(track, itc, 1);
	}
	
	tcbef = track->GetCluster(itcbef);
	tcaft = track->GetCluster(itcaft);
	
	MSG("AlgTrackSR", Msg::kVerbose) << itc
					 << " looking at cluster on plane "
					 << tc->GetPlane() 
					 << " prev cluster on plane "
					 << tcbef->GetPlane()
					 << " netx cluster on plane "
					 << tcaft->GetPlane() << endl;

	//get the slope from the clusters bounding the current one
        localslope = tcaft->GetTPos()-tcbef->GetTPos();
	if(tcaft->GetZPos()!=tcbef->GetZPos())
	  localslope /= tcaft->GetZPos()-tcbef->GetZPos();
// 	else MSG("AlgTrackSR", Msg::kWarning) << "cannot find local slope "
// 					      << "clusters have same zpos  "
// 					      << localslope << endl;
	
	//get the expected tpos at the current cluster
        localtpos = tcbef->GetTPos();
	localtpos += localslope*(tc->GetZPos()-tcbef->GetZPos());

	MSG("AlgTrackSR", Msg::kVerbose) << "local slope = " << localslope
					 << " local tpos = " << localtpos
					 << endl;
	
	
	//find the set of strips in the cluster of size nstripexp which is closest
	//to the expect tpos
	
	//first order the strips by tpos
	for(Int_t i = 0; i<numStrips; ++i){
	  strip = dynamic_cast<CandStripHandle *>(tc->GetStripList()->At(i));
	  stripTPos[i] = strip->GetTPos();
	  stripCharge[i] = strip->GetCharge();
	  MSG("AlgTrackSR", Msg::kVerbose) << stripTPos[i] << endl;
	}

	//sort the array of TPos values in ascending order
	TMath::Sort(numStrips,stripTPos,sortedStrips,false);

        indxbest = 0;
        distbest = 0;
	dist = 0.;
	distsum = 0.;
        for(Int_t indx=0; indx<=numStrips-nstripexp; ++indx){

	  distsum = 0.;

	  for(Int_t istrip=indx; istrip<indx+nstripexp; ++istrip){

	    MSG("AlgTrackSR", Msg::kVerbose) << indx << " " << istrip 
					     << " " << sortedStrips[istrip]
					     << " " << stripTPos[sortedStrips[istrip]]
					     << endl;

	    dist = TMath::Abs(stripTPos[sortedStrips[istrip]]-localtpos);
            if(fParmIsCosmic==1)
              dist *= (0.3+exp(-0.2*stripCharge[sortedStrips[istrip]]));
            
            distsum += dist;
          }
          if (!indx || distsum<distbest) {
            distbest = distsum;
            indxbest = indx;
          }
        }
	
	// count how many hit strips above 1 pe
        nshower=0;
	for(int i= 0; i<=numStrips-nstripexp; ++i){
	  strip = dynamic_cast<CandStripHandle *>(tc->GetStripList()->At(i));
          if(strip->GetCharge()>fParmMinStripPE) ++nshower;
        }

        for(Int_t istrip=indxbest; istrip<indxbest+nstripexp; istrip++) {
          strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(sortedStrips[istrip]));
	  MSG("AlgTrackSR", Msg::kVerbose) << strip->GetPlane() << " " 
					   << strip->GetTPos() << endl;
          cth.AddDaughterLink(*strip);
          cth.SetInShower(strip,nshower-nstripexp);
        }
      }//end if more than expected number of strips in cluster
    }//end loop over track clusters
  }//end loop over tracks


  return;
}

//----------------------------------------------------------------------
Int_t AlgTrackSR::FindNeighborIndex(Track2DSR *track, 
				    Int_t currentIndex,
				    Int_t direction)
{
  MSG("AlgTrackSR", Msg::kVerbose) << "in find neighbor index" << endl;
  
  Bool_t found = false;
  Int_t neighborIndex = currentIndex+direction;
  Int_t nstrip = 0;
  Int_t dstrip = 0;
  TrackClusterSR *tc = 0;

  //loop over clusters in the track starting from the one
  //right next to the current cluster - take the next cluster in the list
  //which has the expected number of strips in it based on the slope at 
  //that cluster
  while(!found 
	&& neighborIndex>=0 && neighborIndex<=track->GetLast()
	){

    MSG("AlgTrackSR", Msg::kVerbose) << "current index = " << currentIndex
				     <<" neighbor index = " << neighborIndex
				     << endl;
    
    tc = track->GetCluster(neighborIndex);
    dstrip = tc->GetStripList()->GetLast()+1;
    dstrip -= TMath::Max(1,(Int_t)(TMath::Abs(track->GetSlope(neighborIndex))/4.));
    if(dstrip<=nstrip) found = true;
          
    ++nstrip;
    neighborIndex += direction;
  }

  neighborIndex -= direction;
  assert(neighborIndex>=0 && neighborIndex<= track->GetLast());
  
  return neighborIndex;
}

//---------------------------------------------------------------------
//figure out if the track is going forward or backwards with respect to
//z.  return the number of points used to determine the direction
Int_t AlgTrackSR::FindTimingDirection(CandTrackSRHandle &cth,
				      CandStripHandleItr &stripItr, 
				      Track2DSR *utrack,
				      Track2DSR *vtrack,
				      Double_t *fitparm,
				      Double_t &maxPathLength,
				      Double_t &timefitchi2)
{

  MSG("AlgTrackSR", Msg::kVerbose) << "in find timing direction" << endl;

  fitparm[0]=0;
  fitparm[1]=0;
  stripItr.Reset();
  Int_t plane = 0;
  maxPathLength = 0.;
  //  Double_t digitpe = 0.;
  // Double_t logadc = 0.;
  const Double_t min_wgt = 0.4;

  CandStripHandle *strip = 0;
  CandDigitHandle *digit = 0;

  StripEnd::StripEnd_t stripEnd = StripEnd::kUnknown;

  //get VldContext and UgliGeomHandle
  CandStripHandle *firststrip = stripItr.Ptr();

  assert(firststrip);

  const VldContext *vldcontext = firststrip->GetVldContext();
  assert(vldcontext);

  UgliGeomHandle ugh(*vldcontext);

  //arrays for determining the timing direction - limit to
  //1000 entries - really shouldnt need more than that to 
  //be able to figure out if the event is going north or 
  //south
  Double_t time[1000];
  Double_t pathLength[1000];
  Double_t weight[1000];
  Int_t digitCtr = 0;

  for(Int_t i = 0; i<1000; ++i){
    time[i] = 0.;
    pathLength[i] = 0.;
    weight[i] = 0.;
  }
  
  //need to put in a check to make sure you dont have more planes in one
  //view than the other

  //get the endpoints for each track
  Int_t uBeg = utrack->GetBegPlane();
  Int_t uEnd = utrack->GetEndPlane();
  Int_t vBeg = vtrack->GetBegPlane();
  Int_t vEnd = vtrack->GetEndPlane();

  Int_t trackBeg = TMath::Min(uBeg,vBeg);
  Int_t trackEnd = TMath::Max(uEnd,vEnd);
  Int_t direction = 1;

  //if the track was loaded up going north to south
  if(uBeg>uEnd){
    direction = -1;
    trackBeg = TMath::Max(uBeg,vBeg);
    trackEnd = TMath::Min(uEnd,vEnd);
  }
  if(TMath::Abs(uBeg-vBeg)>3.){
    if(direction>0) trackBeg = TMath::Max(uBeg,vBeg) - 1;
    else if(direction<0) trackBeg = TMath::Min(uBeg,vBeg) + 1;
  } 
  if(TMath::Abs(uEnd-vEnd)>3.){
    if(direction>0) trackEnd = TMath::Min(uEnd,vEnd) + 1;
    if(direction<0) trackEnd = TMath::Max(uEnd,vEnd) - 1;
  }
  MSG("AlgTrackSR", Msg::kDebug) << "u end points = " << uBeg
				 << " " << uEnd << " v end points = "
				 << vBeg << " " << vEnd << endl
				 << " track end points = " << trackBeg
				 << " " << trackEnd << endl;

  Double_t halfLength = 0.;
  Double_t apos = 0.;
  Double_t distFromCenter = 0.;
  Double_t fiberDist = 0.;
  
  //loop over all strips in the track
  while( (strip = stripItr()) && digitCtr<1000){
    
    plane = strip->GetPlane();

    MSG("AlgTrackSR", Msg::kDebug) << "strip from plane " << plane 
				   << " " << Detector::AsString(fDetector) << endl;

    //only look at double ended strips if in the far detector
    //and strips from planes that have orthogonal measurements 
    if( (strip->GetNDaughters() == 2 || fDetector == Detector::kNear)  
       && plane*direction>=trackBeg*direction 
       && plane*direction<=trackEnd*direction){
      
      MSG("AlgTrackSR", Msg::kVerbose) << "strip good for timing" << endl;

      //get the iterator over the digits in the strip
      CandDigitHandleItr digitItr(strip->GetDaughterIterator());

      while( (digit = digitItr()) && digitCtr<1000){
	
	stripEnd = digit->GetPlexSEIdAltL().GetEnd();
	      
 	pathLength[digitCtr] = cth.GetdS()-cth.GetdS(plane);
	maxPathLength = TMath::Max(maxPathLength, 
				   (double)(cth.GetdS()-cth.GetdS(plane)));
	
	
        //time is recorded in seconds, but i prefer working in ns.
	//get the time from the track rather than the strips
 	time[digitCtr] = cth.GetT(plane,stripEnd)*1.e9;

	//cant use the GetT method for Near detector right now because the near 
	//detector has multiple digits for each strip end and GetT only returns one 
	//time - gotta fix that.
	if(fDetector == Detector::kNear){
	  UgliStripHandle stripHandle = ugh.GetStripHandle(strip->GetStripEndId());
	  halfLength = stripHandle.GetHalfLength();
	  apos = 0.;
	  fiberDist = 0.;
	  
	  if(strip->GetPlaneView() == PlaneView::kV) apos = cth.GetU(plane);
	  else if(strip->GetPlaneView() == PlaneView::kU) apos = cth.GetV(plane);
	  distFromCenter = 0.;
	  if(strip->GetPlaneView() == PlaneView::kV) 
	    distFromCenter = (stripEnd == StripEnd::kNegative) ? apos : -apos;
	  else if(strip->GetPlaneView() == PlaneView::kU) 
	    distFromCenter = (stripEnd == StripEnd::kNegative) ? -apos : apos;
	  
	  fiberDist = (halfLength + distFromCenter + stripHandle.ClearFiber(stripEnd) 
		       + stripHandle.WlsPigtail(stripEnd));
	  
	  time[digitCtr] = strip->GetTime(digit->GetPlexSEIdAltL().GetEnd()) - fiberDist/PropagationVelocity::Velocity();
	  time[digitCtr] *= 1.e9;
	  
	  //digitpe = digit->GetCharge(CalDigitType::kPE);
	  //    if (digitpe>0.0) logadc =  TMath::Log(digitpe)/2.3;
	  //  else {
          //  logadc = 0.0;
	  // MSG("AlgTrackSR", Msg::kWarning)
          //    << "FindTimingDirection() digitpe was "
	  //   << digitpe << ", avoid TMath::Log()." << endl;
	  //  }
	  //	  time[digitCtr] -= (c1*logadc
	  //	     +c2*logadc*logadc
	  //	     +c3*logadc*logadc*logadc);
	}//end if near detector

	//dont use hits in showers for calculating timing
	if(cth.IsInShower(strip)<=fParmMaxClusterSizeTimeFit)
	  weight[digitCtr] = TMath::Min(cth.GetTimeWeight(digit),min_wgt);
	

	MSG("AlgTrackSR", Msg::kDebug) << digitCtr 
				       << " " << pathLength[digitCtr]
				       << " " << time[digitCtr]
				       << " " << strip->GetTime(stripEnd)*1.e9
				       << " " << weight[digitCtr]
				       << endl;
	
	++digitCtr;
	
      }//end loop over digits in strip
    }//end if double ended strip
  }//end loop over strips in track

  Double_t efitparm[2];
  Double_t maxresid = fParmMaxTimeFitRes*1.e9+1.;
  Double_t resid = 0.;
  Int_t ctr = digitCtr;
  Int_t imaxresid = 0;
  Int_t nremove = 0;
  Bool_t dofit=false;
  while(maxresid>fParmMaxTimeFitRes
	&& digitCtr-nremove>fParmMinTimeFitSize 
	&& (1.*nremove)<fParmMaxTimeFitOutlier*digitCtr
	){
    dofit=true;
    LinearFit::Weighted(ctr,pathLength,time,weight,fitparm,efitparm);
    maxresid = 0.;
    imaxresid = 0;
    for(int i=0; i<ctr; ++i){
      resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
      if(weight[i]>0. && TMath::Abs(resid)>maxresid){
        maxresid = TMath::Abs(resid);
        imaxresid = i;
      }
    }
    
    MSG("AlgTrackSR",Msg::kDebug) << "constrained fit, dt/ds slope, offset = " 
				  << fitparm[1] << " " << fitparm[0] << endl;
    MSG("AlgTrackSR",Msg::kDebug) << "maximum residual " << maxresid << " " 
				  << pathLength[imaxresid] << " " 
				  << time[imaxresid] << " " 
				  << weight[imaxresid] << endl;
    
    if(maxresid>fParmMaxTimeFitRes){
      MSG("AlgTrackSR",Msg::kVerbose) << "  removing" << endl;
      weight[imaxresid] = 0.;
      nremove++;
    }
    
    timefitchi2 = 0.;
    for(int i=0; i<ctr; ++i){
      resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
      timefitchi2 += weight[i]*resid*resid;
    }
    MSG("AlgTrackSR",Msg::kDebug) << "chi2 = " << timefitchi2 
				  << "   n = " << digitCtr-nremove << endl;
  }//end loop to find timing fit and remove outliers
  
  timefitchi2 = 0.;
  if(dofit){
    for(int i=0; i<ctr; ++i){
      if(weight[i]>0.){
	MSG("AlgTrackSR",Msg::kDebug) << "TIMEFIT " << i << " " 
				      <<pathLength[i] << " " 
				      << time[i] << " " << weight[i] << endl;
	resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
      timefitchi2 += weight[i]*resid*resid;
      }
    }
  }
    
  MSG("AlgTrackSR",Msg::kDebug) << " TIMEFIT " << fitparm[1] 
				<< " chi^2/ndf = " 
				<< timefitchi2/(1.*(digitCtr-nremove)) 
				<< " " << TMath::Abs(uBeg-uEnd)
				<< " " << TMath::Abs(vBeg-vEnd) << endl;
  
  //check the chi^2 for the fit - if it is really bad dont change the
  //initial guess at the direction for the event, ie just make sure
  //that the slope in time vs pathlength is positive;
  if(fitparm[1]<0. && timefitchi2>=10.*(digitCtr-nremove))
    fitparm[1] *= -1.;

  //return to units of seconds
  fitparm[0] *= 1.e-9;
  fitparm[1] *= 1.e-9;

  return digitCtr-nremove;
}