//------------------------
// String Header First  --
//------------------------
#include <string>

//-----------------------
// This Class's Header --
//-----------------------
#include "TrackingMods/PadTrackMaker.hh"

//---------------
// C++ Headers --
//---------------
#include <algorithm>
#include <iostream>
#include <vector>

//-------------------------------
// Collaborating Class Headers --
//-------------------------------
#include "AbsEnv/AbsEnv.hh"
//#include "HepTuple/HepHist1D.h"
//#include "HepTuple/HepHist2D.h"
//#include "HepTuple/HepNtuple.h"
//#include "HepTuple/HepFileManager.h"
//#include "HepTuple/HepRootFileManager.h"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "TrackingObjects/Storable/CdfTrackColl.hh"
#include "TrackingObjects/SiData/SiHit.hh"
#include "TrackingObjects/Tracks/track_cut.hh"
#include "TrackingObjects/Tracks/CorCotTrack.hh"
#include "VertexObjects/ZVertexColl.hh"
#include "TrackingCT/CT_ReFit.hh"
#include "CotGeometry/CT_Configuration.hh"

//-----------------------------------------------------------------------
// Local Macros, Typedefs, Structures, Unions and Forward Declarations --
//-----------------------------------------------------------------------

static const char rcsid[] = "$Id: PadTrackMaker.cc,v 1.34 2005/08/03 01:28:25 cmircea Exp $";

//----------------
// Constructors --
//----------------
PadTrackMaker::PadTrackMaker( 
    const char* const theName, 
    const char* const theDescription )
    : AppModule( theName, theDescription ),
      _debug("debug", this, false),
      _constrainT0("ConstrainT0",this,true),
      _ptMinForT0Constraint("PTMinForT0Constraint",this,0.3),
      _dropHits("DropHits",this,false),
      _dropHitRoad("DropHitRoad",this,0.06),
      _processName("processName", this, "PROD"),
      _outputCollectionDesc("outputCollectionDesc",this,"PADTracks"),
      _removeSiDuplicates("RemoveSiDuplicates",this,false),
      _cotAxialSegmentsMin("cotAxialSegmentsMin",this,2),
      _cotStereoSegmentsMin("cotStereoSegmentsMin",this,2)
{
  commands()->append(&_debug);
  commands()->append(&_constrainT0);
  commands()->append(&_ptMinForT0Constraint);
  commands()->append(&_dropHits);
  commands()->append(&_dropHitRoad);
  commands()->append(&_processName);
  commands()->append(&_outputCollectionDesc);
  commands()->append(&_removeSiDuplicates); 
  commands()->append(&_cotAxialSegmentsMin);
  commands()->append(&_cotStereoSegmentsMin);
  _debug.addDescription("\tDebug mode");
  _constrainT0.addDescription("\tApply t0 constraint");
  _ptMinForT0Constraint.addDescription("\tMinimum PT for applying t0 constraint");
  _dropHits.addDescription("\tDrop COT hits from fit");
  _dropHitRoad.addDescription("\tRoad for dropping COT hits from fit");
  _processName.addDescription("\tProcess name of view from which PAD track collection will come");
  _removeSiDuplicates.addDescription("\tRemove duplicate IO or SISA tracks");
  _cotAxialSegmentsMin.addDescription("\tMinimum number of COT axial hits");
  _cotStereoSegmentsMin.addDescription("\tMinimum number of COT stereo hits");
  _PionMass = 0.140; 
}

//--------------
// Destructor --
//--------------
PadTrackMaker::~PadTrackMaker( )
{
}

//--------------
// Operations --
//--------------
AppResult PadTrackMaker::beginJob( AbsEvent* anEvent )
{
  // cout << "Hello from Module PadTrackMaker!! Debug value is " <<
  //  _debug.value() << endl;

  _kalFitter = SiKalmanFitter::instance();

  return AppResult::OK;
}

AppResult PadTrackMaker::beginRun( AbsEvent* theRun )
{
  if (_kalFitter){
    _kalFitter->refreshDetector();
  }
  return AppResult::OK;
}

AppResult PadTrackMaker::event( AbsEvent* anEvent )
{
  // iterate over the tracks 
  //
  cotTrkSet.clear();
  CdfTrackView_h hTrackView;
  CdfTrackView::allTracks( hTrackView, _processName.value() );

  ConstHandle<SiHitSet> siliconHitSet_h;

  EventRecord::ConstIterator i( anEvent, "SiHitSet","*", _processName.value());
  for ( ; i.is_valid()&&(i->process_name() == _processName.value()) ; ++i ) 
  {
    siliconHitSet_h = i;
  }
  CdfTrackColl_h padOut = CdfTrackColl_h(new CdfTrackColl(siliconHitSet_h));

  int numtotal=0;
  for (CdfTrackView::const_iterator iTrkLink = hTrackView->contents().begin(); 
       iTrkLink != hTrackView->contents().end(); ++iTrkLink) {
    numtotal++;
  }
  if (_debug.value()) {
    std::cout << "The total number of tracks is " << numtotal << std::endl;
    std::cout << "The Printout of the hTrackView contents follows:" << std::endl;
    std::cout << *hTrackView << std::endl;
  }

  // The real work is done inside makePadTracks here
  if (makePadTracks(anEvent, &*padOut)) {

    if (_debug.value())
      std::cout << "PadTrackMaker -- PAD track collection was made successfully" << std::endl;

    //CdfTrackColl_h padOutH(padOut);
    anEvent->append(padOut);

  }
  // Clear this vector of clnks at the end of the event!
  cotTrkSet.clear();

  if (_debug.value()) {
    std::cout << "The Printout of padOut follows:" << std::endl;
    std::cout << *padOut << std::endl;
    padOut->Print();
  }

  return AppResult::OK;
}

AppResult PadTrackMaker::endJob( AbsEvent* anEvent )
{
  return AppResult::OK;
}

AppModule* PadTrackMaker::clone( const char* cloneName )
{
  return new PadTrackMaker(cloneName,"this module is a clone PadTrackMaker");
}

const char* PadTrackMaker::rcsId( ) const
{
   return rcsid;
}

bool PadTrackMaker::makePadTracks ( EventRecord * anEvent,
				    CdfTrackColl *padOut )
{
  // Create the temporary views (for sorting) and the output collection
  CdfTrackView_h padOIZView  (new CdfTrackView());
  CdfTrackView_h padOISView  (new CdfTrackView());
  CdfTrackView_h padOIView   (new CdfTrackView());
  CdfTrackView_h padCOTView  (new CdfTrackView());
  CdfTrackView_h padSiView   (new CdfTrackView());
  CdfTrackView_h padJunkView (new CdfTrackView());
  CdfTrackView_h padIOView   (new CdfTrackView());

  padOut->set_description(_outputCollectionDesc.value());

  // Loop over the event track collections
  EventRecord::ConstIterator i( anEvent, "CdfTrackColl","*",_processName.value() );
  for ( ; i.is_valid() ; ++i ) 
  {
    if (_debug.value()) {
      std::cout << "FOUND a new collection with process name " << 
      i->process_name() << " and with description " << 
      i->description() << std::endl;
    }    
    if ( i->process_name() != _processName.value() ) 
    {
      if (_debug.value()) {
	std::cout << "  Failed process name check. Continuing..." << std::endl;
      }
      continue;
    }
    if ( i->description() == "GlobalSI_Tracking" || 
	 i->description() == "COT_Global_Tracking" ||
	 i->description() == "IOTracks" )
    {
      ConstHandle<CdfTrackColl> hTrkColl( i );
      const CdfTrackColl::collection_type & trkVect = hTrkColl->contents();
      for ( CdfTrackColl::const_iterator itrk = trkVect.begin(),
	      itrkEnd = trkVect.end();
	    itrk != itrkEnd ; ++itrk )
      {
	//CdfTrack_clnk tempTrack = (*itrk);
	if (_debug.value()) {
	  std::cout << std::endl << "Iterating to a new track inside " <<
	    i->description() << std::endl;
	  std::cout << "It's track " << (*itrk)->id() 
		    << " with algorithm value " 
		    << (*itrk)->algorithm().value()
		    << std::endl;
	}

	// PAD track collection consists only of "leaf" tracks (no children)
	// IO tracks have SI as children. Will keep io tracks and remove their child SI tracks
	if ( !(*itrk)->child() || (*itrk)->algorithm().value()==CdfTrack::InsideOutAlg ) 
	{
	  if (_debug.value()) {
	    std::cout << "Track has no children or is an IO track" << std::endl;
	  }
	  if ( (*itrk)->numCTHits()>0 ) {
	    if ( PassesPadIO((*itrk)) && 
		 (!_removeSiDuplicates.value() || !isDuplicateSi((*itrk),padIOView))) {            
	      if (_debug.value()) {
		std::cout << "Track passed IO selector:"
			  << " pushing onto padIOView" << std::endl;
	      }
	      padIOView->contents().push_back(*itrk);
	    }
	    else if ( PassesPadOIZ(*itrk) ) {
	      if (_debug.value()) {
		std::cout << "Track passed OIZ selector:"
			  << " pushing onto padOIZView." << std::endl;
	      }
	      padOIZView->contents().push_back(*itrk);
	      cotTrkSet.push_back(COTparent(*itrk));
	    }
	    else if ( PassesPadOIS(*itrk) ) { 
	      if (_debug.value()) {
		std::cout << "Track passed OIS selector:"
			  << " pushing onto padOISView." << std::endl;
	      }
	      padOISView->contents().push_back(*itrk);
	      cotTrkSet.push_back(COTparent(*itrk));
	    }
	    // If the two cuts above fail, then we move to the OI
	    // parent track.	    
	    else if ( PassesPadOI(OIparent(*itrk)) ) {
	      if (_debug.value()) {
		std::cout << "Track passed OI selector:"
			  << " pushing onto padOIView" << std::endl;
	      }
	      padOIView->contents().push_back(OIparent(*itrk));
	      cotTrkSet.push_back(COTparent(OIparent(*itrk)));
	    }	    
	    // If the the OI cut fails, then we move to the COT
	    // parent track.
	    else if ( PassesPadCOT(COTparent(*itrk)) ) {
	      if (_debug.value()) {
		std::cout << "Track passed COT selector:"
			  << " pushing onto padCOTView" 
			  << std::endl;
	      }
	      padCOTView->contents().push_back(COTparent(*itrk));
	      cotTrkSet.push_back(COTparent(*itrk));
	    }
	    // The COT tracks which are left go into the junk pile.
	    else {
	      if (_debug.value()) {
		std::cout << "Track couldn't pass any selector: "
			  << " pushing onto padJunkView" << std::endl;
	      }
	      padJunkView->contents().push_back(COTparent(*itrk));
	    }
	  } // has COT hits, is therefore not stand-alone
	  else {
	    // No cuts yet on Si standalone tracks
	    if ( PassesPadSi(*itrk) && 
		 (!_removeSiDuplicates.value() || !isDuplicateSi((*itrk),padSiView))) {
	      if (_debug.value()) {
		std::cout << "Track has no COT hits" << std::endl;
		std::cout << "Track passes Si standalone cut: Pushing onto padSiView" << std::endl;
	      }
	      padSiView->contents().push_back(*itrk);
	    } // passes StandaloneSi
	  } // has COT hits
	} // track has no children 
      } // iterate through tracks
    } // collection passes description cut
  } // loop over collections

  // Now sort the views and copy sequentially the tracks into the 
  // new collection as output.
  std::sort(padOIZView->contents().begin(), padOIZView->contents().end(), DecrPt());
  std::sort(padOISView->contents().begin(), padOISView->contents().end(), DecrPt());
  std::sort(padOIView->contents().begin(), padOIView->contents().end(), DecrPt());
  std::sort(padCOTView->contents().begin(), padCOTView->contents().end(), DecrPt());
  std::sort(padSiView->contents().begin(), padSiView->contents().end(), DecrPt());
  std::sort(padJunkView->contents().begin(), padJunkView->contents().end(), DecrPt());
  std::sort(padIOView->contents().begin(), padIOView->contents().end(), DecrPt());

  // Make copies of the tracks for the new combined output collection.
  // The track algorithms should already be correct.
  // Don't forget the quality flag for the CdfTrack.
  int numParentTracksOut = 0;
  int numOIZTracksOut = 0;
  const CT_Configuration& cfg = CT_Configuration::reference();
  int nsl = cfg.numSuperLayers();
  for (CdfTrackView::const_iterator itrk = padOIZView->contents().begin();
       itrk != padOIZView->contents().end(); ++itrk) {

    CdfTrack* outTrack = new CdfTrack(**itrk);
    outTrack->nullChildLink();
    outTrack->setPADQualityFlag(1);
    
    CdfTrack_clnk oldprecursor = COTparent(outTrack);
    CdfTrack*   precursor    = new CdfTrack(*oldprecursor);
    precursor->setId(oldprecursor->id());

    float constraint0 = -999;    
    if (_constrainT0.value()) {
      // constrain COT parent t0
      EventRecord::ConstIterator iter(anEvent,"ZVertexColl","ZVerCollCotStand",
				      _processName.value());
      
      if (iter.is_valid()) {
	ConstHandle<ZVertexColl> vertexSet(iter);
	ZVertexColl::const_iterator iZVertex;
	for(iZVertex=vertexSet->contents().begin();
	    iZVertex!=vertexSet->contents().end(); iZVertex++) {

	  double vz =iZVertex->zPosition();
	  float t01 =iZVertex->getCotT0();
	  float t02 =iZVertex->getCotT02();
	
	  float trkt0 = precursor->t0();
	  float trkz0 = outTrack->z0();
	  float trkpt = outTrack->pt();

	  float mindz = 99.;
	  float dz = fabs(trkz0 - vz);
	  if (dz < 3. && dz < mindz && trkpt>_ptMinForT0Constraint.value()) {
	    mindz = dz;
	    if (_debug.value()) {
	      std::cout << "z match, track z0: " << 
		trkz0 << ", vertex z: " << vz << std::endl;
	    }
	    if (fabs(trkt0 - t01) < 1.5 && t01 > -99 &&
		fabs(trkt0 - t01) < fabs(trkt0 - t02)) {
	      if (_debug.value()) {
		std::cout << "t0 match, track t0: " << 
		  trkt0 << ", vertex t0: " << t01 << std::endl;
	      }
	      constraint0 = t01;
	    } else if (fabs(trkt0 - t02) < 1.5 && t02 > -99) {
	      if (_debug.value()) {
		std::cout << "t0 match, track t0: " << 
		  trkt0 << ", vertex t0: " << t02 << std::endl;
	      }
	      constraint0 = t02;
	    } else { 
	      if (_debug.value()) {
		std::cout << "no t0 match, track t0: " << 
		  trkt0 << ", vertex t01: " << t01 << 
		  ", vertex t02: " << t02 << std::endl;
	      }
	    }
	  } else {
	    if (_debug.value()) 
	      std::cout << "bad vertex, track z0: " << trkz0 << ", vertex z0: " << vz << std::endl;
	  }
	} // loop over vertices
      } // if valid ZVerCollCotStand
      if (constraint0>-99) {
	CT_ReFit refit;
	refit.useFixed_t0(true);
	refit.setFixedt0(constraint0);
	refit(*precursor);
      
	if (refit.ok()) { 
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    nhits++;
	    if (! ldat.link.is_valid()) {
	      precursor->setResidual(ilr,CT_Residual());
	      continue;
	    }
	    // Save residual
	    precursor->setResidual(ilr,CT_Residual
			  (ldat.hitTV/cfg.driftSpeed(),
			   ldat.hitD,ldat.trackY,ldat.sina));
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setConstrainedT0(constraint0);
	  outTrack->setConstrainedT0(constraint0);
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  if (_debug.value()) std::cout << "refit failed!" << std::endl;
	}
      }
    }

    if (_dropHits.value()) {
      // drop hits with large residuals
      std::vector<int> dropSL;
      std::vector<CT_HitLink> dropHit;
      int ilr = 0;
      for (int isl = 0; isl<nsl; isl++) {
	const CT_Residual* resid = precursor->beginCTResiduals(isl);
	for (int iwr = 0; iwr<cfg.numWires(isl); iwr++) {
	  if (resid[iwr].isValid() && fabs(resid[iwr].deltaD())>_dropHitRoad.value()) {
	    dropSL.push_back(isl);
	    dropHit.push_back(precursor->beginCTHitBlock(isl)[iwr]);
	    precursor->dropHit(ilr);
	    outTrack->dropHit(ilr);
	    if (_debug.value()) 
	      std::cout << "Dropping hit, ilr: " << ilr << ", resid: " << resid[iwr].deltaD() << std::endl;
	  }
	  ilr++;
	}
      }
      if (constraint0>-99) {
	CT_ReFit refit;
	refit.useFixed_t0(true);
	refit.setFixedt0(constraint0);
	refit(*precursor);
      
	if (refit.ok()) { 
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) continue;
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  // restore track
	  for (int nDropped=0; nDropped<dropSL.size(); nDropped++) {
	    precursor->addHit(dropSL[nDropped],dropHit[nDropped]);
	    outTrack->addHit(dropSL[nDropped],dropHit[nDropped]);
	  }
	  if (_debug.value()) {
	    std::cout << "Post hit-drop refit (t0 con) failed!" << std::endl;
	  }
	}
      } else {
	CT_ReFit refit;
	refit(*precursor);
      
	if (refit.ok()) { 
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) continue;
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  // restore track
	  for (int nDropped=0; nDropped<dropSL.size(); nDropped++) {
	    precursor->addHit(dropSL[nDropped],dropHit[nDropped]);
	    outTrack->addHit(dropSL[nDropped],dropHit[nDropped]);
	  }
	  if (_debug.value()) {
	    std::cout << "Post hit-drop refit (no t0 con) failed!" << std::endl;
	  }
	}
      }
    } // drop hits?
	
    CorCotTrack::CorStatus statcor = CorCotTrack::CorTrack(precursor,true);

    const CdfTrack_clnk    oldparent = outTrack->parent();
    if(oldparent.is_nonnull()){
      const CdfTrack_cfwdlnk tempparent(precursor);
      outTrack->setParentLink(tempparent);
    }

    _kalFitter->setCOTErrorScale(1.0);
    _kalFitter->setDEdx(true);
    _kalFitter->setMass(_PionMass);
 
    TrackFitter::FitStatus status = _kalFitter->fit(*outTrack, 0.0);

    if (!(status==TrackFitter::OK)) {
      std::cerr << "Error Kal refitter failed!" << std::endl;
    }
    CorCotTrack::CorStatus uncor = CorCotTrack::UnCorTrack(precursor,true);

    padOut->contents().push_back( outTrack );
    ++numOIZTracksOut;
    padOut->contents().push_back( precursor );
    precursor->setChildLink(outTrack);
    if (_debug.value()) {
      std::cout << "set links: parent: " << precursor->id() << ", child: " << outTrack->id() << std::endl;
    }
    ++numParentTracksOut;
  }

  int numOISTracksOut = 0;
  for (CdfTrackView::const_iterator itrk = padOISView->contents().begin();
       itrk != padOISView->contents().end(); ++itrk) {
    CdfTrack* outTrack = new CdfTrack(**itrk);
    outTrack->nullChildLink();
    drop90ZHits( outTrack );
    outTrack->setPADQualityFlag(1);
    
    CdfTrack_clnk oldprecursor = COTparent(outTrack);
    CdfTrack*   precursor    = new CdfTrack(*oldprecursor);
    precursor->setId(oldprecursor->id());

    float constraint0 = -999;
    if (_constrainT0.value()) {
      // constrain COT parent t0
      EventRecord::ConstIterator iter(anEvent,"ZVertexColl","ZVerCollCotStand",_processName.value());
      if (iter.is_valid()) {
	ConstHandle<ZVertexColl> vertexSet(iter);
	ZVertexColl::const_iterator iZVertex;
	for(iZVertex=vertexSet->contents().begin();
	    iZVertex!=vertexSet->contents().end(); iZVertex++) {
	  
	  double vz =iZVertex->zPosition();
	  float t01 =iZVertex->getCotT0();
	  float t02 =iZVertex->getCotT02();
	
	  float trkt0 = precursor->t0();
	  float trkz0 = outTrack->z0();
	  float trkpt = outTrack->pt();

	  float mindz = 99.;
	  float dz = fabs(trkz0 - vz);
	  if (dz < 3. && dz < mindz && trkpt>_ptMinForT0Constraint.value()) {
	    mindz = dz;
	    if (_debug.value())
	      std::cout << "z match, track z0: " << trkz0 << ", vertex z: " << vz << std::endl;
	    if (fabs(trkt0 - t01) < 1.5 && t01 > -99 &&
		fabs(trkt0 - t01) < fabs(trkt0 - t02)) {
	      if (_debug.value())
		std::cout << "t0 match, track t0: " << trkt0 << ", vertex t0: " << t01 << std::endl;
	      constraint0 = t01;
	    } else if (fabs(trkt0 - t02) < 1.5 && t02 > -99) {
	      if (_debug.value())
		std::cout << "t0 match, track t0: " << trkt0 << ", vertex t0: " << t02 << std::endl;
	      constraint0 = t02;
	    } else {
	      if (_debug.value()) {
		std::cout << "no t0 match, track t0: " << trkt0 << ", vertex t01: " << t01 <<
		  ", vertex t02: " << t02 << std::endl;
	      }
	    }
	  } else {
	    if (_debug.value())
	      std::cout << "bad vertex, track z0: " << trkz0 << ", vertex z0: " << vz << std::endl;
	  }
	} // loop over vertices
      } // if valid ZVerCollCotStand
      if (constraint0>-99) {
	CT_ReFit refit;
	refit.useFixed_t0(true);
	refit.setFixedt0(constraint0);
	refit(*precursor);

	if (refit.ok()) {
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) {
	      precursor->setResidual(ilr,CT_Residual());
	      continue;
	    }
	    // Save residual
	    precursor->setResidual(ilr,CT_Residual
				   (ldat.hitTV/(cfg.driftSpeed()),
				    ldat.hitD,ldat.trackY,ldat.sina));
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setConstrainedT0(constraint0);
	  outTrack->setConstrainedT0(constraint0);
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  if (_debug.value()) std::cout << "refit failed!" << std::endl;
	}
      }
    } // constrain t0?

    if (_dropHits.value()) {
      // drop hits with large residuals
      int ilr = 0;
      std::vector<int> dropSL;
      std::vector<CT_HitLink> dropHit;
      for (int isl = 0; isl<nsl; isl++) {
	const CT_Residual* resid = precursor->beginCTResiduals(isl);
	for (int iwr = 0; iwr<cfg.numWires(isl); iwr++) {
	  if (resid[iwr].isValid() && fabs(resid[iwr].deltaD())>_dropHitRoad.value()) {
	    dropSL.push_back(isl);
	    dropHit.push_back(precursor->beginCTHitBlock(isl)[iwr]);
	    precursor->dropHit(ilr);
	    outTrack->dropHit(ilr);
	    if (_debug.value()) 
	      std::cout << "Dropping hit, ilr: " << ilr << ", resid: " << resid[iwr].deltaD() << std::endl;
	  }
	  ilr++;
	}
      }
      if (constraint0>-99) {
	CT_ReFit refit;
	refit.useFixed_t0(true);
	refit.setFixedt0(constraint0);
	refit(*precursor);
	
	if (refit.ok()) { 
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) continue;
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  // restore track
	  for (int nDropped=0; nDropped<dropSL.size(); nDropped++) {
	    precursor->addHit(dropSL[nDropped],dropHit[nDropped]);
	    outTrack->addHit(dropSL[nDropped],dropHit[nDropped]);
	  }
	  if (_debug.value()) {
	    std::cout << "Post hit-drop refit (t0 con) failed!" << std::endl;
	  }
	}
      } else {
	CT_ReFit refit;
	refit(*precursor);
 
	if (refit.ok()) { 
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) continue;
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  // restore track
	  for (int nDropped=0; nDropped<dropSL.size(); nDropped++) {
	    precursor->addHit(dropSL[nDropped],dropHit[nDropped]);
	    outTrack->addHit(dropSL[nDropped],dropHit[nDropped]);
	  }
	  if (_debug.value()) {
	    std::cout << "Post hit-drop refit (no t0 con) failed!" << std::endl;
	  }
	}
      }
    } // drop hits?

    CorCotTrack::CorStatus statcor = CorCotTrack::CorTrack(precursor,true);

    const CdfTrack_clnk    oldparent = outTrack->parent();
    if(oldparent.is_nonnull()){
      const CdfTrack_cfwdlnk tempparent(precursor);
      outTrack->setParentLink(tempparent);
    }

    _kalFitter->setCOTErrorScale(1.0);
    _kalFitter->setDEdx(true);
    _kalFitter->setMass(_PionMass);
 
    TrackFitter::FitStatus status = _kalFitter->fit(*outTrack, 0.0);

    if (!(status==TrackFitter::OK)) {
      std::cerr << "Error Kal refitter failed!" << std::endl;
    }
    // The refit must be done before the algorithm changes.
    if (outTrack->algorithm().value() == CdfTrack::OutsideIn3DAlg) {
      outTrack->setAlgorithm(CdfTrack::OutsideInStereoAlg);
    }
    else if (outTrack->algorithm().value() == CdfTrack::KalOutsideIn3DAlg) {
      outTrack->setAlgorithm(CdfTrack::KalOutsideInStereoAlg);
    }
    CorCotTrack::CorStatus uncor = CorCotTrack::UnCorTrack(precursor,true);

    padOut->contents().push_back( outTrack );
    ++numOISTracksOut;
    padOut->contents().push_back( precursor );
    precursor->setChildLink(outTrack);
    if (_debug.value()) {
      std::cout << "set links: parent: " << precursor->id() << ", child: " << outTrack->id() << std::endl;
    }
    ++numParentTracksOut;
  }

  int numOITracksOut = 0;
  for (CdfTrackView::const_iterator itrk = padOIView->contents().begin();
       itrk != padOIView->contents().end(); ++itrk) {

    CdfTrack* outTrack = new CdfTrack(**itrk);
    outTrack->nullChildLink();
    drop90ZHits( outTrack );
    dropSASHits( outTrack );
    outTrack->setPADQualityFlag(1);
    
    CdfTrack_clnk oldprecursor = COTparent(outTrack);
    CdfTrack*   precursor    = new CdfTrack(*oldprecursor);
    precursor->setId(oldprecursor->id());

    float constraint0 = -999;
    if (_constrainT0.value()) {
      // constrain COT parent t0
      EventRecord::ConstIterator iter(anEvent,"ZVertexColl","ZVerCollCotStand",_processName.value());
      if (iter.is_valid()) {
	ConstHandle<ZVertexColl> vertexSet(iter);
	ZVertexColl::const_iterator iZVertex;
	for(iZVertex=vertexSet->contents().begin();
	    iZVertex!=vertexSet->contents().end(); iZVertex++) {
	  
	  double vz =iZVertex->zPosition();
	  float t01 =iZVertex->getCotT0();
	  float t02 =iZVertex->getCotT02();

	  float trkt0 = precursor->t0();
	  float trkz0 = outTrack->z0();
	  float trkpt = outTrack->pt();

	  float mindz = 99.;
	  float dz = fabs(trkz0 - vz);
	  if (dz < 6. && dz < mindz && trkpt>_ptMinForT0Constraint.value()) {
	    mindz = dz;
	    if (_debug.value())
	      std::cout << "z match, track z0: " << trkz0 << ", vertex z: " << vz << std::endl;
	    if (fabs(trkt0 - t01) < 1.5 && t01 > -99 &&
		fabs(trkt0 - t01) < fabs(trkt0 - t02)) {
	      if (_debug.value())
		std::cout << "t0 match, track t0: " << trkt0 << ", vertex t0: " << t01 << std::endl;
	      constraint0 = t01;
	    } else if (fabs(trkt0 - t02) < 1.5 && t02 > -99) {
	      if (_debug.value())
		std::cout << "t0 match, track t0: " << trkt0 << ", vertex t0: " << t02 << std::endl;
	      constraint0 = t02;
	    } else {
	      if (_debug.value()) {
		std::cout << "no t0 match, track t0: " << trkt0 << ", vertex t01: " << t01 <<
		  ", vertex t02: " << t02 << std::endl;
	      }
	    }
	  } else {
	    if (_debug.value())
	      std::cout << "bad vertex, track z0: " << trkz0 << ", vertex z0: " << vz << std::endl;
	  }
	} // loop over vertices
      } // if valid ZVerCollCotStand
      if (constraint0>-99) {
	CT_ReFit refit;
	refit.useFixed_t0(true);
	refit.setFixedt0(constraint0);
	refit(*precursor);

	if (refit.ok()) {
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) {
	      precursor->setResidual(ilr,CT_Residual());
	      continue;
	    }
	    // Save residual
	    precursor->setResidual(ilr,CT_Residual
				   (ldat.hitTV/(cfg.driftSpeed()),
				    ldat.hitD,ldat.trackY,ldat.sina));
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setConstrainedT0(constraint0);
	  outTrack->setConstrainedT0(constraint0);
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  if (_debug.value()) std::cout << "refit failed!" << std::endl;
	}
      }
    } // constrain t0?

    if (_dropHits.value()) {
      // drop hits with large residuals
      int ilr = 0;
      std::vector<int> dropSL;
      std::vector<CT_HitLink> dropHit;
      for (int isl = 0; isl<nsl; isl++) {
	const CT_Residual* resid = precursor->beginCTResiduals(isl);
	for (int iwr = 0; iwr<cfg.numWires(isl); iwr++) {
	  if (resid[iwr].isValid() && fabs(resid[iwr].deltaD())>_dropHitRoad.value()) {
	    dropSL.push_back(isl);
	    dropHit.push_back(precursor->beginCTHitBlock(isl)[iwr]);
	    precursor->dropHit(ilr);
	    outTrack->dropHit(ilr);
	    if (_debug.value()) 
	      std::cout << "Dropping hit, ilr: " << ilr << ", resid: " << resid[iwr].deltaD() << std::endl;
	  }
	  ilr++;
	}
      }
      if (constraint0>-99) {
	CT_ReFit refit;
	refit.useFixed_t0(true);
	refit.setFixedt0(constraint0);
	refit(*precursor);
	
	if (refit.ok()) { 
	// compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) continue;
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  // restore track
	  for (int nDropped=0; nDropped<dropSL.size(); nDropped++) {
	    precursor->addHit(dropSL[nDropped],dropHit[nDropped]);
	    outTrack->addHit(dropSL[nDropped],dropHit[nDropped]);
	  }
	  if (_debug.value()) {
	    std::cout << "Post hit-drop refit (t0 con) failed!" << std::endl;
	  }
	}
      } else {
	CT_ReFit refit;
	refit(*precursor);
      
	if (refit.ok()) { 
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) continue;
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  precursor->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  // restore track
	  for (int nDropped=0; nDropped<dropSL.size(); nDropped++) {
	    precursor->addHit(dropSL[nDropped],dropHit[nDropped]);
	    outTrack->addHit(dropSL[nDropped],dropHit[nDropped]);
	  }
	  if (_debug.value()) {
	    std::cout << "Post hit-drop refit (t0 con) failed!" << std::endl;
	  }
	}
      }
    } // drop hits?

    CorCotTrack::CorStatus statcor = CorCotTrack::CorTrack(precursor,true);

    const CdfTrack_clnk    oldparent = outTrack->parent();
    if(oldparent.is_nonnull()){
      const CdfTrack_cfwdlnk tempparent(precursor);
      outTrack->setParentLink(tempparent);
    }

    _kalFitter->setCOTErrorScale(1.0);
    _kalFitter->setDEdx(true);
    _kalFitter->setMass(_PionMass);
 
    TrackFitter::FitStatus status = _kalFitter->fit(*outTrack, 0.0);

    if (!(status==TrackFitter::OK)) {
      std::cerr << "Error Kal refitter failed!" << std::endl;
    }
    if (outTrack->algorithm().value() == CdfTrack::OutsideIn3DAlg) {
      outTrack->setAlgorithm(CdfTrack::OutsideInAlg);
    }
    else if (outTrack->algorithm().value() == CdfTrack::KalOutsideIn3DAlg) {
      outTrack->setAlgorithm(CdfTrack::KalOutsideInAlg);
    }
    CorCotTrack::CorStatus uncor = CorCotTrack::UnCorTrack(precursor,true);

    padOut->contents().push_back( outTrack );
    ++numOITracksOut;
    padOut->contents().push_back( precursor );
    precursor->setChildLink(outTrack);
    ++numParentTracksOut;
  }

  int numIOTracksOut = 0;
  for (CdfTrackView::const_iterator itrk = padIOView->contents().begin();
       itrk != padIOView->contents().end(); ++itrk) {
    if (isDuplicateIO(*itrk)) continue;
    CdfTrack* outTrack = new CdfTrack(**itrk);
    outTrack->nullChildLink();
    outTrack->setPADQualityFlag(1);  
    padOut->contents().push_back( outTrack );
    ++numIOTracksOut;
  }

  int numCOTTracksOut = 0;
  for (CdfTrackView::const_iterator itrk = padCOTView->contents().begin();
       itrk != padCOTView->contents().end(); ++itrk) {
    CdfTrack* outTrack = new CdfTrack(**itrk);
    outTrack->nullChildLink();
    outTrack->setPADQualityFlag(1);

    float constraint0 = -999;
    if (_constrainT0.value()) {
      // constrain COT parent t0
      EventRecord::ConstIterator iter(anEvent,"ZVertexColl","ZVerCollCotStand",
				      _processName.value());
      if (iter.is_valid()) {
	ConstHandle<ZVertexColl> vertexSet(iter);
	ZVertexColl::const_iterator iZVertex;
	for(iZVertex=vertexSet->contents().begin();
	    iZVertex!=vertexSet->contents().end(); iZVertex++) {

	  double vz =iZVertex->zPosition();
	  float t01 =iZVertex->getCotT0();
	  float t02 =iZVertex->getCotT02();

	  float trkt0 = outTrack->t0();
	  float trkz0 = outTrack->z0();
	  float trkpt = outTrack->pt();
	
	  float mindz = 99.;
	  float dz = fabs(trkz0 - vz);
	  if (dz < 6. && dz < mindz && trkpt>_ptMinForT0Constraint.value()) {
	    mindz = dz;
	    if (_debug.value())
	      std::cout << "z match, track z0: " << trkz0 << ", vertex z: " << vz << std::endl;
	    if (fabs(trkt0 - t01) < 1.5 && t01 > -99 &&
		fabs(trkt0 - t01) < fabs(trkt0 - t02)) {
	      if (_debug.value())
		std::cout << "t0 match, track t0: " << trkt0 << ", vertex t0: " << t01 << std::endl;
	      constraint0 = t01;
	    } else if (fabs(trkt0 - t02) < 1.5 && t02 > -99) {
	      if (_debug.value())
		std::cout << "t0 match, track t0: " << trkt0 << ", vertex t0: " << t02 << std::endl;
	      constraint0 = t02;
	    } else {
	      if (_debug.value()) {
		std::cout << "no t0 match, track t0: " << trkt0 << ", vertex t01: " << t01 <<
		  ", vertex t02: " << t02 << std::endl;
	      }
	    }
	  } else {
	    if (_debug.value())
	      std::cout << "bad vertex, track z0: " << trkz0 << ", vertex z0: " << vz << std::endl;
	  }
	} // loop over vertices
      } // if valid ZVerCollCotStand
      if (constraint0>-99) {
	CT_ReFit refit;
	refit.useFixed_t0(true);
	refit.setFixedt0(constraint0);
	refit(*outTrack);

	if (refit.ok()) {
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) {
	      outTrack->setResidual(ilr,CT_Residual());
	      continue;
	    }
	    // Save residual
	    outTrack->setResidual(ilr,CT_Residual
				  (ldat.hitTV/(cfg.driftSpeed()),
				   ldat.hitD,ldat.trackY,ldat.sina));
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  outTrack->setConstrainedT0(constraint0);
	  outTrack->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  if (_debug.value()) std::cout << "refit failed!" << std::endl;
	}
      }
    } // constrain t0?
    
    if (_dropHits.value()) {
      // drop hits with large residuals
      int ilr = 0;
      std::vector<int> dropSL;
      std::vector<CT_HitLink> dropHit;
      for (int isl = 0; isl<nsl; isl++) {
	const CT_Residual* resid = outTrack->beginCTResiduals(isl);
	for (int iwr = 0; iwr<cfg.numWires(isl); iwr++) {
	  if (resid[iwr].isValid() && fabs(resid[iwr].deltaD())>_dropHitRoad.value()) {
	    dropSL.push_back(isl);
	    dropHit.push_back(outTrack->beginCTHitBlock(isl)[iwr]);
	    outTrack->dropHit(ilr);
	    if (_debug.value()) 
	      std::cout << "Dropping hit, ilr: " << ilr << ", resid: " << resid[iwr].deltaD() << std::endl;
	  }
	  ilr++;
	}
      }
      if (constraint0>-99) {
	CT_ReFit refit;
	refit.useFixed_t0(true);
	refit.setFixedt0(constraint0);
	refit(*outTrack);
      
	if (refit.ok()) { 
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) continue;
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  outTrack->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  // restore track
	  for (int nDropped=0; nDropped<dropSL.size(); nDropped++) 
	    outTrack->addHit(dropSL[nDropped],dropHit[nDropped]);
	  
	  if (_debug.value()) {
	    std::cout << "Post hit-drop refit (t0 con) failed!" << std::endl;
	  }
	}
      } else {
	CT_ReFit refit;
	refit(*outTrack);
	
	if (refit.ok()) { 
	  // compute chi2
	  const int nlr = cfg.layer0(nsl - 1) + cfg.numWires(nsl - 1);
	  double chi2 = 0.;
	  int nhits = 0;
	  for (int ilr = 0; ilr < nlr; ++ilr) {
	    const CT_ReFit::LayerData& ldat = *(refit.layer_data() + ilr);
	    if (! ldat.link.is_valid()) continue;
	    nhits++;
	    chi2 += (ldat.deltaD()*ldat.deltaD()) * ldat.weight;
	  }
	  outTrack->setFit3D(refit.npar(), refit.par(), refit.cov(),chi2);
	} else {
	  // restore track
	  for (int nDropped=0; nDropped<dropSL.size(); nDropped++) 
	    outTrack->addHit(dropSL[nDropped],dropHit[nDropped]);

	  if (_debug.value()) {
	    std::cout << "Post hit-drop refit (t0 con) failed!" << std::endl;
	  }
	}
      }
    } // drop hits?

    padOut->contents().push_back( outTrack );
    ++numCOTTracksOut;
  }

  int numSiTracksOut = 0;
  for (CdfTrackView::const_iterator itrk = padSiView->contents().begin();
       itrk != padSiView->contents().end(); ++itrk) {
    CdfTrack* outTrack = new CdfTrack(**itrk);
    outTrack->nullChildLink();
    outTrack->setPADQualityFlag(1);
    padOut->contents().push_back( outTrack );
    ++numSiTracksOut;
  }

  int numJunkTracksOut = 0;
   // Need to set a quality flag here to distinguish
  for (CdfTrackView::const_iterator itrk = padJunkView->contents().begin();
       itrk != padJunkView->contents().end(); ++itrk) {
    CdfTrack* outTrack = new CdfTrack(**itrk);
    outTrack->nullChildLink();
    padOut->contents().push_back( outTrack );
    ++numJunkTracksOut;
  }

  // the last part of the collection is for IO parents
  for (CdfTrackView::const_iterator itrk = padIOView->contents().begin();
       itrk != padIOView->contents().end(); ++itrk) {
    if (isDuplicateIO(*itrk)) continue;
    CdfTrack_clnk theIOcot = COTparent(*itrk);
    CdfTrack* ultimateParent = new CdfTrack( *theIOcot );
    padOut->contents().push_back( ultimateParent );
    ++numParentTracksOut;
  }
    
  if (_debug.value()) {
    std::cout << "There were " << numOIZTracksOut << " OIZ tracks written out."
	      << std::endl;
    std::cout << "There were " << numOISTracksOut << " OIS tracks written out."
	      << std::endl;
    std::cout << "There were " << numOITracksOut << " OI tracks written out."
	      << std::endl;
    std::cout << "There were " << numIOTracksOut << " IO tracks written out."
	      << std::endl;
    std::cout << "There were " << numCOTTracksOut << " COT tracks written out."
	      << std::endl;
    std::cout << "There were " << numSiTracksOut << " Si tracks written out."
	      << std::endl;
    std::cout << "There were " << numJunkTracksOut << " junk tracks written out."
	      << std::endl;
    std::cout << "There were " << numParentTracksOut << " OI parent tracks written out." << std::endl;
  }

  // Assign track ID's in the new collection so the links can be made
  padOut->assignTrackIds( anEvent );

  // Set parent-child links for the silicon tracks and their COT parents
  /* need to think about this */
  int numTracks = 0;
  int numOITracks = numOIZTracksOut + numOISTracksOut + numOITracksOut ;
  int numTracksBeforeIOParents = 2*numOITracks + numIOTracksOut +
    numCOTTracksOut + numSiTracksOut + numJunkTracksOut;

  CdfTrackColl::const_iterator firstParentTrk, firstIOTrk;
  CdfTrackColl::const_iterator iend = padOut->contents().end();
  for ( CdfTrackColl::iterator itrk = padOut->contents().begin();
	itrk != iend; ++itrk) {

    // find start of IO
    if (numTracks == 2*numOITracks) {
      firstIOTrk = itrk;
    }
    if (numTracks == numTracksBeforeIOParents) {
      firstParentTrk = itrk;
    }
    ++numTracks;
  }

  // first do OIZ, OIS, OI
  int numSiTracksLinked = 0;
  int numSiliTracksOut =  numOITracks; 
  CdfTrackColl::const_iterator ptrk = padOut->contents().begin() + 1;
  for ( CdfTrackColl::const_iterator itrk = padOut->contents().begin();
	itrk != iend && numSiTracksLinked < numSiliTracksOut; 
	itrk=itrk+2, ptrk=ptrk+2) {
    (*itrk)->setParentLink(CdfTrack_clnk(ConstLink<CdfTrackColl>(padOut),ptrk));
    (*ptrk)->setChildLink (CdfTrack_clnk(ConstLink<CdfTrackColl>(padOut),itrk));
    ++numSiTracksLinked;
  }

  // then do IO
  int numSiliconTracksLinked = 0;
  int numSiliconTracksOut =  numIOTracksOut;
  CdfTrackColl::const_iterator jtrk = firstParentTrk;
  for ( CdfTrackColl::const_iterator itrk = firstIOTrk;
        itrk != iend && numSiliconTracksLinked < numSiliconTracksOut;
        ++itrk, ++jtrk) {
    (*itrk)->setParentLink(CdfTrack_clnk(ConstLink<CdfTrackColl>(padOut),jtrk));
    (*jtrk)->setChildLink (CdfTrack_clnk(ConstLink<CdfTrackColl>(padOut),itrk));
    ++numSiliconTracksLinked;
  }

  return true;
} // makePadTracks
  

CdfTrack_clnk PadTrackMaker::COTparent(CdfTrack_clnk inputTrack) {
  CdfTrack_clnk parentTrack = inputTrack;
  while (parentTrack->algorithm().value() !=
	 CdfTrack::CotStandAloneAlg && parentTrack->parent()) { 
    // This loop ends when parentTrack is a COT track or if
    // parentTrack is the ultimate parent.
    if (_debug.value()) {
      std::cout << "COTparent: current track alg id is " <<
	parentTrack->algorithm().value() << std::endl;
    }
    parentTrack = parentTrack->parent();
  }
  if (parentTrack->algorithm().value() == CdfTrack::CotStandAloneAlg) {
    return parentTrack;
  }
  else {
    if (_debug.value()) {
      std::cout << "Error!  Here's a COT parent that's not COT!" 
		<< "  Returning input track." 
		<< std::endl;
    }
    return inputTrack;
  }
}

CdfTrack_clnk PadTrackMaker::OIparent(CdfTrack_clnk inputTrack) {
  CdfTrack_clnk parentTrack = inputTrack;
  while (parentTrack->algorithm().value() != CdfTrack::OutsideInAlg &&
	 parentTrack->algorithm().value() != CdfTrack::KalOutsideIn3DAlg &&
	 parentTrack->parent()) { 
    // This loop ends when parentTrack is a OI track or if
    // parentTrack is the ultimate parent.
    if (_debug.value()) {
      std::cout << "OIparent: current track alg id is " <<
	parentTrack->algorithm().value() << std::endl;
    }
    parentTrack = parentTrack->parent();
  }
  if (parentTrack->algorithm().value() == CdfTrack::OutsideInAlg ||
      parentTrack->algorithm().value() == CdfTrack::KalOutsideIn3DAlg ) {
    return parentTrack;
  }
  else {
    if (_debug.value()) {
      std::cout << "Error!  Here's a OI parent that's not OI!" 
		<< "  Returning input track." 
		<< std::endl;
    }
    return inputTrack;
  }
}

bool PadTrackMaker::PassesPadOIZ( CdfTrack_clnk trk ) {
  // Don't throw away a great silicon track just because the COT hits
  // are lacking...

  // first find number of good COT superlayers
  int nst=0;
  int nax=0;
  // stereo first
  for (int sl=0; sl<8; sl=sl+2) 
    if (trk->numCTHits(sl)>=5) nst++;
  // now axial
  for (int sl=1; sl<8; sl=sl+2)
    if (trk->numCTHits(sl)>=5) nax++;

  if ( (trk->algorithm().value() == CdfTrack::OutsideIn3DAlg ||
	trk->algorithm().value() == CdfTrack::KalOutsideIn3DAlg) &&
       ( (
	  ((nax>= _cotAxialSegmentsMin.value()&& nst>=_cotStereoSegmentsMin.value()) || 
	   (nax>=_cotAxialSegmentsMin.value() && nst>=1 && trk->lastLayerCT()<95 )) &&
	  trk->numSIHitsPhi()>=3 && trk->numSIHitsPhi()<=5 ) || 
	 trk->numSIHitsPhi()>5 ) &&
       ( (trk->numSIHitsSt()==0 && trk->numSIHitsZ()>=3) ||
	 (trk->numSIHitsSt()>0  && trk->numSIHitsZ()>0 ) ) &&
       sqrt(trk->getHelixFit().getErrorMatrix()(3,3))<0.05 ) {
    return true;
  }
  else {
    return false;
  }
}

bool PadTrackMaker::PassesPadOIS( CdfTrack_clnk trk ) {

  // first find number of good COT superlayers
  int nst=0;
  int nax=0;
  // stereo first
  for (int sl=0; sl<8; sl=sl+2) 
    if (trk->numCTHits(sl)>=5) nst++;
  // now axial
  for (int sl=1; sl<8; sl=sl+2)
    if (trk->numCTHits(sl)>=5) nax++;

  if ( (trk->algorithm().value() == CdfTrack::OutsideIn3DAlg ||
	trk->algorithm().value() == CdfTrack::KalOutsideIn3DAlg) &&
       ( (
	  ((nax>=_cotAxialSegmentsMin.value() && nst>=_cotStereoSegmentsMin.value()) || 
	   (nax>=_cotAxialSegmentsMin.value() && nst>=1 && trk->lastLayerCT()<95 )) &&
	  trk->numSIHitsPhi()>=3 && trk->numSIHitsPhi()<=5 ) || 
	 trk->numSIHitsPhi()>5 ) &&
       trk->numSIHitsSt()>=1 &&
       sqrt(trk->getHelixFit().getErrorMatrix()(3,3))<0.5 ) {
    return true;
  }
  else {
    return false;
  }
}

bool PadTrackMaker::PassesPadOI( CdfTrack_clnk trk ) {

  // first find number of good COT superlayers
  int nst=0;
  int nax=0;
  // stereo first
  for (int sl=0; sl<8; sl=sl+2) 
    if (trk->numCTHits(sl)>=5) nst++;
  // now axial
  for (int sl=1; sl<8; sl=sl+2)
    if (trk->numCTHits(sl)>=5) nax++;

  if ( (trk->algorithm().value() == CdfTrack::OutsideInAlg ||
	trk->algorithm().value() == CdfTrack::KalOutsideIn3DAlg) &&
       ((
	 ((nax>=_cotAxialSegmentsMin.value() && nst>=_cotStereoSegmentsMin.value()) || 
	  (nax>=_cotAxialSegmentsMin.value() && nst>=1 && trk->lastLayerCT()<95 )) &&
	 trk->numSIHitsPhi()>=3 && trk->numSIHitsPhi()<=5 ) || 
	trk->numSIHitsPhi()>5 ) &&
       sqrt(trk->getHelixFit().getErrorMatrix()(4,4))<0.05 ) {
    return true;
  }
  else {
    return false;
  }
}

bool PadTrackMaker::PassesPadCOT( CdfTrack_clnk trk ) {
  // first find number of good COT superlayers
  int nst=0;
  int nax=0;
  // stereo first
  for (int sl=0; sl<8; sl=sl+2) 
    if (trk->numCTHits(sl)>=5) nst++;
  // now axial
  for (int sl=1; sl<8; sl=sl+2)
    if (trk->numCTHits(sl)>=5) nax++;

  if ( trk->algorithm().value() == CdfTrack::CotStandAloneAlg &&
       ((nax>=_cotAxialSegmentsMin.value() && nst>=_cotStereoSegmentsMin.value()) || 
	(nax>=_cotAxialSegmentsMin.value() && nst>=1 && trk->lastLayerCT()<95 )) &&
       trk->chi2CT()/trk->numCTHits() < 10.) {
    return true;
  }
  else {
    return false;
  }
}

bool PadTrackMaker::PassesPadSi( CdfTrack_clnk trk ) {
  if ( (trk->algorithm().value() == CdfTrack::KalSvxStandAloneAlg || trk->algorithm().value() == CdfTrack::KalForwardStandAloneAlg) && 
      trk->numSIHitsPhi()>=3) {
 
    if (_debug.value() && trk->parent().is_nonnull())
      std::cout << "SI parent algorithm is " << trk->parent()->algorithm().value() << std::endl;

    if (trk->lastLayerCT()<(CT_Configuration::reference()).numLayers()-1 &&
	trk->parent().is_null())
      return true;
  }
 
  return false;
}

bool PadTrackMaker::PassesPadIO(CdfTrack_clnk trk) {

  // first find number of good COT superlayers
  int nst=0;
  int nax=0;

  if (trk->algorithm().value() == CdfTrack::InsideOutAlg && 
      trk->lastLayerCT() == 95) {
  // stereo first
    for (int sl=0; sl<8; sl=sl+2) 
      if (trk->numCTHits(sl)>=5) nst++;
    // now axial
    for (int sl=1; sl<8; sl=sl+2)
      if (trk->numCTHits(sl)>=5) nax++;
  }

  if ( trk->algorithm().value() == CdfTrack::InsideOutAlg &&
       ((nax>=_cotAxialSegmentsMin.value() && nst>=_cotStereoSegmentsMin.value()) || 
	(trk->lastLayerCT()<95 )))
    return true;
  else 
    return false;
}

bool PadTrackMaker::isDuplicateIO(CdfTrack_clnk trk) {
  // Check for COT overlap
  const CdfTrack_clnk par = trk->parent();
  if (verbose()) printf("check if io is a duplicate. cotTrkSet size=%d\n",cotTrkSet.size());
  for (TrackSet::iterator i = cotTrkSet.begin(); i!=cotTrkSet.end(); ++i) {
    float f = TrackSetComp::commonHitFraction(*i,par);
    if (verbose()) printf("%.3f\t", f);
    if (f>TrackSetComp::GetOverlap()) return true;
  }
  if (verbose()) printf("\n");
  return false;
}

bool PadTrackMaker::isDuplicateSi(CdfTrack_clnk trk, CdfTrackView_h trackView) {
  // Check for Silicon Overlap
  
  for (CdfTrackView::const_iterator itrk = trackView->contents().begin();
       itrk != trackView->contents().end(); ++itrk) {
    if (trk->usedSI()== (*itrk)->usedSI()) {
      bool matchedHit = false;
      bool allMatchedHits = true;
      for(CdfTrack::SiHitIterator hitIter = trk->beginSIHits();  
	  hitIter != trk->endSIHits(); ++hitIter){
	matchedHit =  false;
	for(CdfTrack::SiHitIterator hitIter2 = (*itrk)->beginSIHits();  
	    hitIter2 != (*itrk)->endSIHits(); ++hitIter2){
	  if ((*hitIter)->id() == (*hitIter2)->id())   matchedHit = true;
	}
	if (!matchedHit) allMatchedHits = false;
      }
      if (allMatchedHits){ 
	if (verbose()){
	  std::cout << "Found SiSA duplicate: " << std::endl;
	  trk->print();
	  (*itrk)->print();
	}
	return true;
      }
    }
  }
  return false;
}

void PadTrackMaker::drop90ZHits( CdfTrack* trk ) {
  std::vector<const SiHit*> hitsToKill;
  for(CdfTrack::SiHitIterator hitIter = trk->beginSIHits();
      hitIter != trk->endSIHits(); ++hitIter) {
    if ( (*hitIter)->nSideHit() ) {
      int layer = (*hitIter)->getSiDigiCode().getLayer();
      if ( layer==1 || layer==2 || layer==4) {
	hitsToKill.push_back(*hitIter);
      }
    }
  }
  for (int i=0; i<hitsToKill.size(); i++) {
    trk->dropHit(hitsToKill[i]);
  }
}

void PadTrackMaker::dropSASHits( CdfTrack* trk) {
  std::vector<const SiHit*> hitsToKill;
  for(CdfTrack::SiHitIterator hitIter = trk->beginSIHits();
      hitIter != trk->endSIHits(); ++hitIter) {
    if ( (*hitIter)->nSideHit() ) {
      int layer = (*hitIter)->getSiDigiCode().getLayer();
      if ( layer==3 || layer==5 || layer==6 || layer==7 ) {
	hitsToKill.push_back(*hitIter);
      }
    }
  }
  for (int i=0; i<hitsToKill.size(); i++) {
    trk->dropHit(hitsToKill[i]);
  }
}