DTrackCandidate_factory.cc

Go to the documentation of this file.

// $Id$
//
//    File: DTrackCandidate_factory.cc
// Created: Mon Jul 18 15:23:04 EDT 2005
// Creator: davidl (on Darwin wire129.jlab.org 7.8.0 powerpc)
//

/// Factory to match track candidates generated by the FDC- and CDC-specific 
/// track finding codes.  The code uses several algorithms in the following 
/// order of precedence:
///      Method 1:  Swims cdc candidates pointing in the forward direction 
///                 and tries to match them geometrically to fdc candidates.
///      Method 2:  Projects the fdc candidates into the cdc region and tries
///                 to match to the cdc hits in each cdc candidate that points 
///                 in the forward direction.
///      Method 3:  Similar to method 2, except that it projects a cdc 
///                 candidate into the FDC region.
///      Method 4:  Matches left-over fdc candidates to other left-over fdc 
///                 candidates.
///      Method 5:  Matches cdc candidates that do not point toward the cdc
///                 endplate that were not matched by previous methods to fdc
///                 candidates.
///      Method 6:  Tries to match stray unused cdc hits with fdc candidates
///      Method 7:  Alternate method to match leftover fdc candidates that 
///                 have already been matched to other track candidates
///      Method 8:  Attempts to "improve" a cdc candidate to be matched to an
///                 FDC candidate with an assumption of the hit position in 
///                 outermost stereo straw
///      Methods 9-13:  Various tricks to try to associate segments in the FDC
///                 with other segments or tracks
/// If a match is found, the code attempts to improve the track parameters by
/// redoing the Riemann Circle fit with the additional hits.  If an fdc 
/// candidate is matched to a cdc candidate, or several previously unused hits
/// are added to the track, the code attempts to place the "start" position 
/// parameters of the track at a radius just outside the start counter.

#include <cmath>
using namespace std;

#include "DTrackCandidate_factory.h"
#include "DTrackCandidate_factory_CDC.h"
#include "START_COUNTER/DSCHit.h"
#include "DANA/DApplication.h"
#include <JANA/JCalibration.h>
#include <TRACKING/DHoughFind.h>

#include <TROOT.h>
#include <TH2F.h>
#include <TMath.h>

#define CUT 10.
#define RADIUS_CUT 50.0
#define BEAM_VAR 1. // cm^2
#define EPS 0.001

//------------------
// cdc_fdc_match
//------------------
inline bool cdc_fdc_match(double p_fdc,double p_cdc,double dist){
  //double frac=fabs(1.-p_cdc/p_fdc);
  //double frac2=fabs(1.-p_fdc/p_cdc);
  double p=p_fdc;
  if (p_cdc <p ) p=p_cdc;
  if (dist<10. && dist <4.+1.75/p
      //&& (frac<0.5 || frac2<0.5)
      ) return true;
  return false;
}

//------------------
// SegmentSortByLayerincreasing
//------------------
inline bool SegmentSortByLayerincreasing(const DFDCSegment* const &segment1, const DFDCSegment* const &segment2) {
   // Compare DFDCSegment->DFDCPseudo[0]->DFDCWire->layer
   int layer1 = 100; // defaults just in case there is a segment with no hits
   int layer2 = 100;
   
   if(segment1->hits.size()>0)layer1=segment1->hits[0]->wire->layer;
   if(segment2->hits.size()>0)layer2=segment2->hits[0]->wire->layer;

   return layer1 < layer2;
}

//------------------
// CDCHitSortByLayerincreasing
//------------------
inline bool CDCHitSortByLayerincreasing(const DCDCTrackHit* const &hit1, const DCDCTrackHit* const &hit2) {
   // Used to sort CDC hits by layer (ring) with innermost layer hits first

   // if same ring, sort by wire number
   if(hit1->wire->ring == hit2->wire->ring)
   {
      if(hit1->wire->straw == hit2->wire->straw)
         return (hit1->dE > hit2->dE);
      return hit1->wire->straw < hit2->wire->straw;
   }

   return hit1->wire->ring < hit2->wire->ring;
}

//------------------
// FDCHitSortByLayerincreasing
//------------------
inline bool FDCHitSortByLayerincreasing(const DFDCPseudo* const &hit1, const DFDCPseudo* const &hit2) {
  // Used to sort FDC hits by layer with most upstream layer hits first  
  // if same layer, sort by wire number
  if(hit1->wire->layer == hit2->wire->layer)
    {
      if(hit1->wire->wire == hit2->wire->wire)
   return (hit1->dE > hit2->dE);
      return (hit1->wire->wire < hit2->wire->wire);
    }
  
  return hit1->wire->layer < hit2->wire->layer;
}

//------------------
// init
//------------------
jerror_t DTrackCandidate_factory::init(void)
{
   MAX_NUM_TRACK_CANDIDATES = 20;

   return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DTrackCandidate_factory::brun(JEventLoop* eventLoop,int32_t runnumber){
  DApplication* dapp=dynamic_cast<DApplication*>(eventLoop->GetJApplication());
  const DGeometry *dgeom  = dapp->GetDGeometry(runnumber);
  bfield = dapp->GetBfield(runnumber);
  FactorForSenseOfRotation=(bfield->GetBz(0.,0.,65.)>0.)?-1.:1.;

  // Get start counter geometry;
  dgeom->GetStartCounterGeom(sc_pos,sc_norm);

  dIsNoFieldFlag = (dynamic_cast<const DMagneticFieldMapNoField*>(bfield) != NULL);
  if(dIsNoFieldFlag)
  {
    //Setting this flag makes it so that JANA does not delete the objects in _data.  This factory will manage this memory.
      //This is all of these pointers are just copied from the "StraightLine" factory, and should not be re-deleted.
    SetFactoryFlag(NOT_OBJECT_OWNER);
  }
  else
    ClearFactoryFlag(NOT_OBJECT_OWNER); //This factory will create it's own objects.

  // Get the position of the exit of the CDC endplate from DGeometry
  double endplate_z,endplate_dz,endplate_rmin;
  dgeom->GetCDCEndplate(endplate_z,endplate_dz,endplate_rmin,endplate_rmax);
  cdc_endplate.SetZ(endplate_z+endplate_dz);

  dParticleID = NULL;
  eventLoop->GetSingle(dParticleID);

  JCalibration *jcalib = dapp->GetJCalibration(runnumber);
  map<string, double> targetparms;
  if (jcalib->Get("TARGET/target_parms",targetparms)==false){
    TARGET_Z = targetparms["TARGET_Z_POSITION"];
   }
  else{
    dgeom->GetTargetZ(TARGET_Z);
  }

   // Initialize the stepper
  stepper=new DMagneticFieldStepper(bfield);
  stepper->SetStepSize(1.0);

  DEBUG_HISTS=false;
  gPARMS->SetDefaultParameter("TRKFIND:DEBUG_HISTS",DEBUG_HISTS);

  if(DEBUG_HISTS){
    dapp->Lock();
    /*
    match_center_dist2=(TH2F*)gROOT->FindObject("match_center_dist2");
    if (!match_center_dist2){
      match_center_dist2=new TH2F("match_center_dist2","larger radius vs matching distance squared between two circle centers",100,0,100.,100,0,100);
      match_center_dist2->SetYTitle("r_{c} [cm]");
      match_center_dist2->SetXTitle("(#Deltad)^{2} [cm^{2}]");
    }
    */
    match_dist=(TH2F*)gROOT->FindObject("match_dist");
    if (!match_dist){
      match_dist=new TH2F("match_dist","Matching distance",
           120,0.,60.,500,0,25.);
      match_dist->SetXTitle("r (cm)");
      match_dist->SetYTitle("#Deltar (cm)");
    }
    match_dist_vs_p=(TH2F*)gROOT->FindObject("match_dist_vs_p");
    if (!match_dist_vs_p) {
      match_dist_vs_p=new TH2F("match_dist_vs_p","Matching distance vs p",
                50,0.,7.,100,0,25.);
      match_dist_vs_p->SetYTitle("#Deltar (cm)");
      match_dist_vs_p->SetXTitle("p (GeV/c)");
    }
    dapp->Unlock();
  }

  gPARMS->SetDefaultParameter("TRKFIND:MAX_NUM_TRACK_CANDIDATES", MAX_NUM_TRACK_CANDIDATES); 

  //  MIN_NUM_HITS=6;
  //gPARMS->SetDefaultParameter("TRKFIND:MIN_NUM_HITS", MIN_NUM_HITS);
  
  DEBUG_LEVEL=0;
  gPARMS->SetDefaultParameter("TRKFIND:DEBUG_LEVEL", DEBUG_LEVEL);

  return NOERROR;
}

//------------------
// erun
//------------------
jerror_t DTrackCandidate_factory::erun(void)
{
  if (stepper) {
    delete stepper;
    stepper = nullptr;
  }

  return NOERROR;
}

//------------------
// erun
//------------------
jerror_t DTrackCandidate_factory::fini(void)
{
  if (stepper) {
    delete stepper;
    stepper = nullptr;
  }

  return NOERROR;
}



//------------------
// evnt
//------------------
jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
{
  if(dIsNoFieldFlag)
  {
    //Clear previous objects: //JANA doesn't do it because NOT_OBJECT_OWNER was set
     //It DID delete them though, in the "StraightLine" factory
    _data.clear();

    vector<const DTrackCandidate*> locStraightLineCandidates;
    loop->Get(locStraightLineCandidates, "StraightLine");
    for(size_t loc_i = 0; loc_i < locStraightLineCandidates.size(); ++loc_i)
      _data.push_back(const_cast<DTrackCandidate*>(locStraightLineCandidates[loc_i]));
    return NOERROR;
  }
  
   // Start counter hits
  vector<const DSCHit *>schits;
  loop->Get(schits);


  // Clear private vectors
  cdctrackcandidates.clear();
  fdctrackcandidates.clear();
  trackcandidates.clear();
  mycdchits.clear();

  // Get the track candidates from the CDC and FDC candidate finders
  loop->Get(cdctrackcandidates, "CDC");
  loop->Get(fdctrackcandidates, "FDCCathodes");

  // List of cdc hits
  loop->Get(mycdchits);

  // Vectors for keeping track of cdc matches and projections to the endplate
  vector<unsigned int> cdc_forward_ids;
  vector<unsigned int> cdc_backward_ids;
  vector<DVector3> cdc_endplate_projections;
 
  // Vector to keep track of cdc hits used in candidates
  vector<unsigned int>used_cdc_hits(mycdchits.size());

  // vector to keep track of the matching status of each fdc candidate
  vector<int>forward_matches(fdctrackcandidates.size());
   
  // Normal vector for CDC endplate
  DVector3 norm(0,0,1);

  // Keep track of CDC hits that have already been used in track candidates
  for(unsigned int i=0; i<cdctrackcandidates.size(); i++){ 
    const DTrackCandidate *srccan = cdctrackcandidates[i];
    for (unsigned int n=0;n<srccan->used_cdc_indexes.size();n++){
      used_cdc_hits[srccan->used_cdc_indexes[n]]=1;
    }
  }
  unsigned int num_unmatched_cdcs=0;
  for (unsigned int i=0;i<used_cdc_hits.size();i++){
    if (used_cdc_hits[i]==0) num_unmatched_cdcs++;
  }

  //Loop through the list of CDC candidates, flagging those that point toward 
  // the FDC.
  for(unsigned int i=0; i<cdctrackcandidates.size(); i++){  
    const DTrackCandidate *srccan = cdctrackcandidates[i];
    DVector3 mom=srccan->momentum();
    DVector3 pos=srccan->position();
    
    // Propagate track to CDC endplate
    bool isForward=false;
    if (fdctrackcandidates.size()>0){
      if (mom.Theta()<M_PI_2){
   // First do a quick projection using a helical model to see if it is 
   // worth adding this cdc candidate to the list of forward-going tracks
   // that could pass into the FDC...
   ProjectHelixToZ(cdc_endplate.z(),srccan->charge(),mom,pos);
   
   if (pos.Perp()<60.){
     // do an actual swim to the cdc endplate
     mom=srccan->momentum();
     pos=srccan->position();
     stepper->SetCharge(srccan->charge());
     stepper->SwimToPlane(pos,mom,cdc_endplate,norm,NULL);
     if (pos.Perp()<60.){
       cdc_endplate_projections.push_back(pos);
       cdc_forward_ids.push_back(i);
       isForward=true;
     }
   }
      }
    }
    if (isForward==false){
      cdc_backward_ids.push_back(i);
    }
  }

  // Variables for candidate number accounting
  int num_forward_cdc_cands_remaining=cdc_forward_ids.size();
  int num_fdc_cands_remaining=fdctrackcandidates.size();
  vector<int>cdc_forward_matches(cdc_forward_ids.size()); 
  vector<int>cdc_backward_matches(cdc_backward_ids.size());

  // Loop through the list of FDC candidates looking for matches between the
  // CDC and the FDC in the transition region.
  if (num_forward_cdc_cands_remaining>0){
    for(unsigned int i=0; i<fdctrackcandidates.size(); i++){
      // Break out if there are no forward-pointing cdc candidates that have 
      // not already been matched to fdc candidates.
      if (num_forward_cdc_cands_remaining==0) break;
      
      // The current FDC candidate
      const DTrackCandidate *fdccan = fdctrackcandidates[i];
      vector<const DFDCSegment *>segments;
      fdccan->GetT(segments);
      stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);
      if (segments[0]->package!=0) continue;

      // Flag for matching status
      bool got_match=false;
    
      if (DEBUG_LEVEL>0) _DBG_ << "Attempting FDC/CDC matching method #1..." <<endl; 

      // Check that the z-position at the doca to the beam line is within the 
      // active volume of the detector to prevent connecting low momentum 
      // curlers coming off the cdc endplate at a small angle with CDC
      // candidates
      if (CheckZPosition(fdccan)==false){
   // mark to avoid trying to match to CDC with the alternate 
   // matching algorithms below
   forward_matches[i]=-1;
   continue;
      }

      // First try the matching method that projects the cdc candidate to the
      // cdc endplate where the fdc candidates are reported.
      if (MatchMethod1(fdccan,cdc_forward_ids,cdc_endplate_projections,
             cdc_forward_matches)){
   if (DEBUG_LEVEL>0) _DBG_ << "... matched to FDC candidate #" << i <<endl;

   got_match=true;
      }
      else{  
   if (DEBUG_LEVEL>0) _DBG_ << "Attempting FDC/CDC matching method #2..." <<endl; 
   // loop over the forward-pointing cdc candidates
   for (unsigned int j=0;j<cdc_forward_ids.size();j++){
     if (cdc_forward_matches[j]) continue;
     const DTrackCandidate *cdccan=cdctrackcandidates[cdc_forward_ids[j]];
     
     // Try to gather up stray CDC hits from candidates that were not 
     // matched with the previous algorithm by projecting the helical track
     // from the fdc into the cdc region
     if (MatchMethod2(fdccan,cdccan)){
       if (DEBUG_LEVEL>0) _DBG_ << "... matched to FDC candidate #" << i <<endl;
       // mark the cdc candidate as matched
       cdc_forward_matches[j]=1;
       
       if (DEBUG_LEVEL>0)
         _DBG_ << ".. matched to CDC candidate #" <<cdc_forward_ids[j] <<endl;
       got_match=true;
       break;
     }
   }
      }
  
      if (got_match){
   // Mark the FDC candidate as matched
   forward_matches[i]=1;
   num_fdc_cands_remaining--;
   
   // update number of unmatched forward-pointing cdc candidates
   num_forward_cdc_cands_remaining--;  
      }
    } // loop over fdc candidates
  } // check for forward-pointing cdc candidates
  
  // If starting with the fdc candidates did not lead to a complete set of
  // CDC-FDC matches, try looping over the remaining CDC candidates that point
  // toward the FDC.
  if (num_forward_cdc_cands_remaining>0){
    for (unsigned int i=0;i<cdc_forward_ids.size();i++){ 
      if (cdc_forward_matches[i]) continue;
      if (num_forward_cdc_cands_remaining==0) break;

      if (DEBUG_LEVEL>0) _DBG_ << "Attempting FDC/CDC matching method #3..." <<endl; 

      // This method projects the cdc track into the FDC region
      const DTrackCandidate *cdccan=cdctrackcandidates[cdc_forward_ids[i]];
      if (MatchMethod3(cdccan,forward_matches)){
   num_fdc_cands_remaining--;
   num_forward_cdc_cands_remaining--;
   
   //Mark the cdc track candidate as matched
   cdc_forward_matches[i]=1;

   if (DEBUG_LEVEL>0) _DBG_ << "... matched to CDC candidate #" << i <<endl;
      }
    } // loop over forward-pointing cdc candidates
  }

  // The following uses a trick to improve the circle fit and "fix" the 
  // direction of a cdc candidate for those candidates whose outermost hit 
  // is a stereo straw by assuming that the z position is at the end of this
  // straw at the cdc endplate, thereby giving us an additional point in the
  // xy plane that can be used in the circle fit.  The code then attempts to 
  // match this "improved" cdc candidate with the remaining fdc candidates.
  if (num_fdc_cands_remaining>0 && num_forward_cdc_cands_remaining>0){
    for (unsigned int j=0;j<cdc_forward_ids.size();j++){
      if (num_fdc_cands_remaining==0) break;
      if (num_forward_cdc_cands_remaining==0) break;
      if (cdc_forward_matches[j]==0){ 
   const DTrackCandidate *cdccan = cdctrackcandidates[cdc_forward_ids[j]];
   
        if (MatchMethod8(cdccan,forward_matches)==true){
     cdc_forward_matches[j]=1;
     num_fdc_cands_remaining--;
     num_forward_cdc_cands_remaining--;
   }

      } 
    } 
  }

  // add to the main list of candidates those we did not successfully merge
  if (num_forward_cdc_cands_remaining>0){
    for (unsigned int j=0;j<cdc_forward_ids.size();j++){
      if (cdc_forward_matches[j]==0){
   DTrackCandidate *can = new DTrackCandidate;

   const DTrackCandidate *cdccan = cdctrackcandidates[cdc_forward_ids[j]];
   vector<const DCDCTrackHit *>cdchits;
   cdccan->GetT(cdchits);
   
   // circle parameters
   can->rc=cdccan->rc;
   can->xc=cdccan->xc;
   can->yc=cdccan->yc;
   
   DVector3 mom=cdccan->momentum();
   DVector3 pos=cdccan->position();
   
   can->setMomentum(mom);
   can->setPosition(pos);
   can->setPID(cdccan->PID());
   
   for (unsigned int n=0;n<cdchits.size();n++){
     used_cdc_hits[cdccan->used_cdc_indexes[n]]=1;
     can->AddAssociatedObject(cdchits[n]);
   }
   can->chisq=cdccan->chisq;
   can->Ndof=cdccan->Ndof;
   trackcandidates.push_back(can);
      }  
    }
  }

  for (unsigned int j=0;j<cdc_backward_ids.size();j++){    
    const DTrackCandidate *cdccan = cdctrackcandidates[cdc_backward_ids[j]]; 

    // Sometimes the cdc track candidate parameters for tracks that are actually
    // going forward are so poor that they don't seem to point towards the cdc
    // end plate, so the previous matching methods fail.  Try one more time 
    // to match these to FDC segments by using Method 8...
    if (num_fdc_cands_remaining>0 
   && MatchMethod8(cdccan,forward_matches)==true){
      num_fdc_cands_remaining--;
      cdc_backward_matches[j]=1;
    }
    else{
      DVector3 mom=cdccan->momentum();
      DVector3 pos=cdccan->position();

      // Check for candidates that appear to go backwards but are actually 
      // going forwards and try to match these to remaining fdc candidates
      if (num_fdc_cands_remaining>0 && mom.Theta()>M_PI_2 && !sc_pos.empty()){
   if (TryToFlipDirection(schits,mom,pos)){
     if (DEBUG_LEVEL>0){
       _DBG_<< "Flipped direction of track (backward to forward) to look for FDC match..." << endl;
     }
     DVector3 norm(0.,0.,1.);
     double p_cdc=mom.Mag();
     stepper->SetCharge(cdccan->charge());
     stepper->SwimToPlane(pos,mom,cdc_endplate,norm,NULL);

     for (unsigned int i=0;i<fdctrackcandidates.size();i++){
       if (num_fdc_cands_remaining==0) break;
       if (forward_matches[i]==0){
         const DTrackCandidate *fdccan=fdctrackcandidates[i];
         if (MatchMethod2(fdccan,cdccan)){
      if (DEBUG_LEVEL>0) {
        _DBG_ << "... matched to FDC candidate #" << i <<endl;
      }
      forward_matches[i]=1;
      cdc_backward_matches[j]=1;
      num_fdc_cands_remaining--;
      break;
         }

         // Use MatchMethod1 if the previous method did not work
         DVector3 fdcpos=fdccan->position();
         DVector3 fdcmom=fdccan->momentum();
         double p_fdc=fdcmom.Mag();
         ProjectHelixToZ(cdc_endplate.z(),fdccan->charge(),fdcmom,fdcpos);
         double diff=(pos-fdcpos).Mag(); 
         if (cdc_fdc_match(p_fdc,p_cdc,diff)){
      // Get the segment data
      vector<const DFDCSegment *>segments;
      fdccan->GetT(segments);
      
      // The following code snippet is intended to prevent matching a cdc
      // track candidate to a low momentum particle curling from the cdc endplate  
      double theta=fdcmom.Theta();
      if (segments.size()>1 && p_fdc<0.3 && theta< 5.*M_PI/180.){ 
        continue;
      }
      
      if (MakeCandidateFromMethod1(theta,segments,cdccan)){
        forward_matches[i]=1;
        num_fdc_cands_remaining--;
        cdc_backward_matches[j]=1;
        if (DEBUG_LEVEL>0)
          _DBG_ << ".. matched to CDC candidate #" << cdc_backward_ids[j] <<endl;
        break;
      }
         } // match criterion met?
       } // fdc candidate not already matched?
     } // loop over fdc candidates
   }
      } // do we have fdc candidates remaining?
    }
  } //loop over "backward" going cdc tracks

  // put the remaining "backward" tracks in the main list of candidates
  for (unsigned int j=0;j<cdc_backward_ids.size();j++){
    if (cdc_backward_matches[j]==0){
      const DTrackCandidate *cdccan = cdctrackcandidates[cdc_backward_ids[j]]; 

      // Get the cdc hits
      vector<const DCDCTrackHit *>cdchits;
      cdccan->GetT(cdchits);
      stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing);
     
      DTrackCandidate *can = new DTrackCandidate;
      can->setMomentum(cdccan->momentum());
      can->setPosition(cdccan->position());
      can->setPID(cdccan->PID());
      
      // circle parameters
      can->rc=cdccan->rc;
      can->xc=cdccan->xc;
      can->yc=cdccan->yc;
      
      for (unsigned int n=0;n<cdchits.size();n++){
   can->AddAssociatedObject(cdchits[n]);
      }
      can->chisq=cdccan->chisq;
      can->Ndof=cdccan->Ndof;
      
      trackcandidates.push_back(can);
    }
  }

  // If there are more than one FDC candidates remaining, use the best track 
  // candidate knowledge we have so far to recheck for matches to other fdc 
  // candidates.
  if (num_fdc_cands_remaining>1){ 
    for (unsigned int j=0;j<fdctrackcandidates.size();j++){
      if (num_fdc_cands_remaining<1) break;
      if (forward_matches[j]<=0){
   // Try to match to fdc candidates that have not already been matched
   // to another track
   if (MatchMethod4(fdctrackcandidates[j],forward_matches,num_fdc_cands_remaining)==true){
     forward_matches[j]=1;
     num_fdc_cands_remaining--;
   }
      }
    }
  }

  // Of the remaining fdc candidates, output those that have more than one 
  // segment associated with them.
  if (num_fdc_cands_remaining>0){ 
    for (unsigned int j=0;j<fdctrackcandidates.size();j++){
      if (forward_matches[j]<=0){
   const DTrackCandidate *fdccan = fdctrackcandidates[j];
   
   // Get the segment data
   vector<const DFDCSegment *>segments;
   fdccan->GetT(segments);
   
   if (segments.size()>1){
     // Create new candidate for output
     DTrackCandidate *can = new DTrackCandidate;
     // circle parameters
     can->rc=fdccan->rc;
     can->xc=fdccan->xc;
     can->yc=fdccan->yc;
     
     can->setMomentum(fdccan->momentum());
     can->setPosition(fdccan->position());
     can->setPID(fdccan->PID());
     can->chisq=fdccan->chisq;
     can->Ndof=fdccan->Ndof; 
     for (unsigned int m=0;m<segments.size();m++){      
       for (unsigned int n=0;n<segments[m]->hits.size();n++){
         const DFDCPseudo *fdchit=segments[m]->hits[n];
         can->AddAssociatedObject(fdchit);
       }
     }
     trackcandidates.push_back(can);
     
     num_fdc_cands_remaining--;
     forward_matches[j]=1;
   }
      }
    }
  }

  // Try to match remaining fdc candidates to track candidates that have 
  // track segments that have already been matched together
  if (num_fdc_cands_remaining>0){
    vector<int>candidate_updated(trackcandidates.size());
    for (unsigned int i=0;i<trackcandidates.size();i++){
      if (num_fdc_cands_remaining<=0) break;
      
      DTrackCandidate *can=trackcandidates[i];      
      if (MatchMethod7(can,forward_matches,num_fdc_cands_remaining)==false){
   if (can->momentum().Mag()<0.5){
     if (MatchMethod12(can,forward_matches,num_fdc_cands_remaining)){
       candidate_updated[i]=1;
     }
   }
      }
    }
    if (num_fdc_cands_remaining>0){
      // if the candidate was updated, try again...
      for (unsigned int i=0;i<trackcandidates.size();i++){
   if (num_fdc_cands_remaining==0) break;
   if (candidate_updated[i]==1){
     DTrackCandidate *can=trackcandidates[i];   
     if (MatchMethod7(can,forward_matches,num_fdc_cands_remaining)==false){
       if (can->momentum().Mag()<0.5){
         MatchMethod12(can,forward_matches,num_fdc_cands_remaining);
       }
     }
   }
      }
    }
  }

  // We should be left with only single-segment fdc candidates.  We try to 
  // connect them together using alternate fitting techniques, either by
  // forcing the circle to originate from the target or at the other extreme
  // not including the fake point at the origin.
  if (num_fdc_cands_remaining){
    bool matched_stray_segments=MatchStraySegments(forward_matches,
                     num_fdc_cands_remaining);
    if (num_fdc_cands_remaining && matched_stray_segments){
      for (unsigned int i=0;i<trackcandidates.size();i++){
   if (num_fdc_cands_remaining==0) break;
   
   DTrackCandidate *can=trackcandidates[i];      
   if (MatchMethod7(can,forward_matches,num_fdc_cands_remaining)==false){
     if (can->momentum().Mag()<0.5){
       MatchMethod12(can,forward_matches,num_fdc_cands_remaining);
     }
   }
      } 
    }
    if (num_fdc_cands_remaining){ 
      // Not much we can do here -- add to the final list of candidates
      for (unsigned int j=0;j<forward_matches.size();j++){
   if (num_fdc_cands_remaining==0) break;
   if (forward_matches[j]<=0){
     const DTrackCandidate *srccan=fdctrackcandidates[j];
     // Get the segment data
     vector<const DFDCSegment *>segments;
     srccan->GetT(segments);
     const DFDCSegment *segment=segments[0];

     DTrackCandidate *can = new DTrackCandidate;
     // circle parameters
     can->rc=srccan->rc;
     can->xc=srccan->xc;
     can->yc=srccan->yc;

     can->Ndof=srccan->Ndof;
     can->chisq=srccan->chisq;
     can->setPID(srccan->PID());
     can->setMomentum(srccan->momentum());
     can->setPosition(srccan->position());
     for (unsigned int n=0;n<segment->hits.size();n++){
       const DFDCPseudo *fdchit=segment->hits[n];
       can->AddAssociatedObject(fdchit);
     }
     
     trackcandidates.push_back(can);
   }
      }
    }
  }

  //There are frequently several hits in the CDC that are not associated with
  // any track segment.  In this case, try to attach these stray hits to a
  // track in the FDC that has hits in the first package.
  if (num_unmatched_cdcs>0){
    if (DEBUG_LEVEL>0){
      _DBG_ << "Trying to use stray CDC hits in match to FDC..." << endl;
    }

    for (unsigned int i=0;i<trackcandidates.size();i++){
      if (num_unmatched_cdcs>0){
   DTrackCandidate *candidate=trackcandidates[i];
   // Get the hits from the candidate
       vector<const DCDCTrackHit *>usedcdchits;
       candidate->GetT(usedcdchits);
       if (usedcdchits.size()==0){
    vector<const DFDCPseudo*>usedfdchits;
    candidate->GetT(usedfdchits);
    // Sort the hits 
    stable_sort(usedfdchits.begin(),usedfdchits.end(),FDCHitSortByLayerincreasing);
    if ((usedfdchits[0]->wire->layer-1)/6==0){
      MatchMethod6(candidate,usedfdchits,used_cdc_hits,num_unmatched_cdcs);
    }
       }
      }
    }
  } 

  // Only output the candidates that have at least a minimum number of hits
  for (unsigned int i=0;i<trackcandidates.size();i++){
    DTrackCandidate *candidate=trackcandidates[i];
     // Get the hits from the candidate
    vector<const DFDCPseudo*>myfdchits;
    candidate->GetT(myfdchits);
    vector<const DCDCTrackHit *>mycdchits;
    candidate->GetT(mycdchits);
    
    if (mycdchits.size()>=3 || myfdchits.size()>=3){
      _data.push_back(candidate);
    }
    else delete candidate;
  }

  if((int(_data.size()) > MAX_NUM_TRACK_CANDIDATES) && (MAX_NUM_TRACK_CANDIDATES >= 0))
   {
     if (DEBUG_LEVEL>0) _DBG_ << "Number of candidates = " << _data.size()
               << " > " <<  MAX_NUM_TRACK_CANDIDATES
               << " --- skipping track fitting! " << endl;
     for(size_t loc_i = 0; loc_i < _data.size(); ++loc_i)
       delete _data[loc_i];
     _data.clear();
   }


  // Set CDC ring & FDC plane hit patterns
  for(size_t loc_i = 0; loc_i < _data.size(); ++loc_i)
  {
    vector<const DCDCTrackHit*> locCDCTrackHits;
    _data[loc_i]->Get(locCDCTrackHits);

    vector<const DFDCPseudo*> locFDCPseudos;
    _data[loc_i]->Get(locFDCPseudos);

    _data[loc_i]->dCDCRings = dParticleID->Get_CDCRingBitPattern(locCDCTrackHits);
    _data[loc_i]->dFDCPlanes = dParticleID->Get_FDCPlaneBitPattern(locFDCPseudos);
  }
  
  return NOERROR;
}


// Obtain position and momentum at the exit of a given package using the 
// helical track model.
void DTrackCandidate_factory::GetPositionAndMomentum(const DFDCSegment *segment,
                       DVector3 &pos, 
                       DVector3 &mom) const{
  // Position of track segment at last hit plane of package
  double x=segment->xc+segment->rc*cos(segment->Phi1);
  double y=segment->yc+segment->rc*sin(segment->Phi1);
  double z=segment->hits[0]->wire->origin.z();

  // Track parameters
  //double kappa=segment->q/(2.*segment->rc);
  double phi0=segment->phi0;
  double tanl=segment->tanl;
  double z0=segment->z_vertex;

  // Useful intermediate variables
  double cosp=cos(phi0);
  double sinp=sin(phi0);
  double sperp=(z-z0)/tanl;
  //double twoks=2.*kappa*sperp;
  double twoks=FactorForSenseOfRotation*segment->q*sperp/segment->rc;
  double sin2ks=sin(twoks);
  double cos2ks=cos(twoks); 

  // Get Bfield
  double B=fabs(bfield->GetBz(x,y,z));

  // Momentum
  double pt=0.003*B*segment->rc;
  double px=pt*(cosp*cos2ks-sinp*sin2ks);
  double py=pt*(sinp*cos2ks+cosp*sin2ks);
  double pz=pt*tanl;

  pos.SetXYZ(x,y,z);
  mom.SetXYZ(px,py,pz);
}

// Get position and momentum at doca to beam line
void DTrackCandidate_factory::GetPositionAndMomentum(const DHelicalFit &fit,
                       double Bz,
                       DVector3 &pos,
                       DVector3 &mom) const{
  // Find position at doca to beam line
  double phi0=atan2(-fit.x0,fit.y0);
  if (fit.h<0) phi0+=M_PI;
  double sinphi0=sin(phi0);
  double sign=(sinphi0>0)?1.:-1.;
  if (fabs(sinphi0)<1e-8) sinphi0=sign*1e-8;
  double cosphi0=cos(phi0);
  double D=FactorForSenseOfRotation*fit.h*fit.r0-fit.x0/sinphi0;
  double x=-D*sinphi0;
  double y=D*cosphi0;
  double dx=pos.x()-x;
  double dy=pos.y()-y;
  double ratio=sqrt(dx*dx+dy*dy)/(2.*fit.r0);
  double phi_s=(ratio<1.)?2.*asin(ratio):M_PI;
  double newz=pos.z()-phi_s*fit.tanl*fit.r0;
  pos.SetXYZ(x,y,newz);
    
  // momentum at POCA to beam line
  double pt=0.003*Bz*fit.r0;
  mom.SetXYZ(pt*cosphi0,pt*sinphi0,pt*fit.tanl);
}

// Get the position and momentum at a fixed radius from the beam line
jerror_t DTrackCandidate_factory::GetPositionAndMomentum(DHelicalFit &fit,
                      double Bz,
                      const DVector3 &origin,
                      DVector3 &pos,
                      DVector3 &mom) const{
  double r2=90.0;
  double xc=fit.x0;
  double yc=fit.y0;
  double rc=fit.r0;
  double tworc=2.*rc;
  double rc2=rc*rc;
  double xc2=xc*xc;
  double yc2=yc*yc;
  double xc2_plus_yc2=xc2+yc2;
  double a=(r2-xc2_plus_yc2-rc2)/tworc;
  double b=xc2_plus_yc2-a*a;
  if (b<0){
    // We did not find an intersection between the two circles, so return 
    // an error.  The values of mom and pos are not changed. 
    return VALUE_OUT_OF_RANGE;
  }

  double temp1=yc*sqrt(b);
  double temp2=xc*a;
  double cosphi_plus=(temp2+temp1)/xc2_plus_yc2;
  double cosphi_minus=(temp2-temp1)/xc2_plus_yc2;

  // Direction tangent and transverse momentum
  double tanl=fit.tanl;
  double pt=0.003*Bz*rc;

  double phi_plus=acos(cosphi_plus);
  double phi_minus=acos(cosphi_minus);
  double x_plus=xc+rc*cosphi_plus;
  double x_minus=xc+rc*cosphi_minus;
  double y_plus=yc+rc*sin(phi_plus);
  double y_minus=yc+rc*sin(phi_minus);

  // if the resulting radial position on the circle from the fit does not agree
  // with the radius to which we are matching, we have the wrong sign for phi+ 
  // or phi-
  double r2_plus=x_plus*x_plus+y_plus*y_plus;
  double r2_minus=x_minus*x_minus+y_minus*y_minus;  
  if (fabs(r2-r2_plus)>EPS){
    phi_plus*=-1.;
    y_plus=yc+rc*sin(phi_plus);
  }
  if (fabs(r2-r2_minus)>EPS){
    phi_minus*=-1.;
    y_minus=yc+rc*sin(phi_minus);
  }

  // Choose phi- or phi+ depending on proximity to one of the cdc hits
  double xwire=origin.x();
  double ywire=origin.y();
  double dx=x_minus-xwire;
  double dy=y_minus-ywire;
  double d2_minus=dx*dx+dy*dy;
  dx=x_plus-xwire;
  dy=y_plus-ywire;
  double d2_plus=dx*dx+dy*dy;
 
  DVector3 pos0(pos); // save the input position, for use in finding z
  if (d2_plus>d2_minus){
    fit.h=-1.;
    phi_minus=M_PI-phi_minus;
    pos.SetXYZ(x_minus,y_minus,0.); // z will be filled later
    mom.SetXYZ(pt*sin(phi_minus),pt*cos(phi_minus),pt*tanl);
  }
  else{
    fit.h=1.;
    phi_plus*=-1.;   
    pos.SetXYZ(x_plus,y_plus,0.); // z will be filled later
    mom.SetXYZ(pt*sin(phi_plus),pt*cos(phi_plus),pt*tanl);
  }
  // Next find the z-position corresponding to the new (x,y) position
  double ratio=(pos0-pos).Perp()/tworc;
  double sperp=(ratio<1.)?tworc*asin(ratio):tworc*M_PI_2;
  pos.SetZ(pos0.z()-sperp*tanl);
  
  return NOERROR;
}

// Find the position along a helical path at the z-position z
void DTrackCandidate_factory::ProjectHelixToZ(const double z,const double q,
                     const DVector3 &mom,
                     DVector3 &pos){
  double pt=mom.Perp();
  double phi=mom.Phi();
  double sinphi=sin(phi);
  double cosphi=cos(phi);
  double tanl=tan(M_PI_2-mom.Theta());
  double x0=pos.X();
  double y0=pos.Y();
  double z0=pos.Z();
  double B=fabs(bfield->GetBz(x0,y0,z0));
  double twokappa=FactorForSenseOfRotation*0.003*B*q/pt;
  double one_over_twokappa=1./twokappa;
  double sperp=(z-z0)/tanl;
  double twoks=twokappa*sperp;
  double sin2ks=sin(twoks);
  double one_minus_cos2ks=1.-cos(twoks);
  double x=x0+(cosphi*sin2ks-sinphi*one_minus_cos2ks)*one_over_twokappa;
  double y=y0+(sinphi*sin2ks+cosphi*one_minus_cos2ks)*one_over_twokappa;
  
  pos.SetXYZ(x,y,z);
}


// Routine to do a crude match between cdc wires and a helical approximation to
// the trajectory
double DTrackCandidate_factory::DocaToHelix(const DCDCTrackHit *hit,
                   double q,
                   const DVector3 &pos,
                   const DVector3 &mom){
  double pt=mom.Perp();
  double phi=mom.Phi();
  double sinphi=sin(phi);
  double cosphi=cos(phi);
  double tanl=tan(M_PI_2-mom.Theta());
  double x0=pos.X();
  double y0=pos.Y();
  double z0=pos.Z();
  //  _DBG_ <<endl;
  double B=fabs(bfield->GetBz(x0,y0,z0));
  double twokappa=FactorForSenseOfRotation*0.003*B*q/pt;
  double one_over_twokappa=1./twokappa;
  
  double sperp=0;
  DVector3 origin=hit->wire->origin;
  double z0w=origin.z();
  DVector3 dir=(1./hit->wire->udir.z())*hit->wire->udir;
  
  DVector3 wirepos;
  double old_doca2=1e8;
  double doca2=old_doca2;
  double z=z0;
  double x=x0,y=y0;
  while (z>50. && z<180. && x<60. && y<60.){
    old_doca2=doca2;
    sperp-=1.;
    double twoks=twokappa*sperp;
    double sin2ks=sin(twoks);
    double one_minus_cos2ks=1.-cos(twoks);
    x=x0+(cosphi*sin2ks-sinphi*one_minus_cos2ks)*one_over_twokappa;
    y=y0+(sinphi*sin2ks+cosphi*one_minus_cos2ks)*one_over_twokappa;
    z=z0+sperp*tanl;
    
    wirepos=origin+(z-z0w)*dir;
    double dxw=x-wirepos.x();
    double dyw=y-wirepos.y();
    
    doca2=dxw*dxw+dyw*dyw;
    if (doca2>old_doca2){
      break;
    } 
  }

  return old_doca2;
}


// Find the particle charge given the helical fit result and an fdc hit 
// on the track
double DTrackCandidate_factory::GetSenseOfRotation(DHelicalFit &fit,
                 const DFDCPseudo *fdchit,
                 const DVector3 &pos
                 ){
  // Get circle parameters
  double rc=fit.r0;
  double xc=fit.x0;
  double yc=fit.y0;
  // Compute phi rotation from "vertex" to fdc hit
  double dphi=(fdchit->wire->origin.Z()-pos.z())/(rc*fit.tanl);
  double Phi1=atan2(pos.Y()-yc,pos.X()-xc);
  
  // Positive and negative changes in phi
  double phiplus=Phi1+dphi;
  double phiminus=Phi1-dphi;
  DVector2 plus(xc+rc*cos(phiplus),yc+rc*sin(phiplus));
  DVector2 minus(xc+rc*cos(phiminus),yc+rc*sin(phiminus));

  // Compute differences 
  double d2plus=(plus-fdchit->xy).Mod2();
  double d2minus=(minus-fdchit->xy).Mod2();

  if (d2minus<d2plus)
    return (-1.0);
  else
    return (+1.0);
}



// Redo the circle fit using all of the cdc axial wires and fdc hits associated
// with the track candidate.  Also compute the average Bz.
jerror_t DTrackCandidate_factory::DoRefit(DHelicalFit &fit,
                 vector<const DFDCSegment *>segments,
                 vector<const DCDCTrackHit *>cdchits,
                 double &Bz){
  unsigned int num_hits=0;
  // Initialize Bz
  Bz=0.;

  // Add the cdc axial wires to the list of hits to use in the fit 
  for (unsigned int k=0;k<cdchits.size();k++){  
    if (cdchits[k]->is_stereo==false){
      double cov=0.213;  //guess
      const DVector3 origin=cdchits[k]->wire->origin;
      double x=origin.x(),y=origin.y(),z=origin.z();
      fit.AddHitXYZ(x,y,z,cov,cov,0.,true);
      Bz+=bfield->GetBz(x,y,z);
      num_hits++;
    }   
  }
  // Add the FDC hits and estimate Bz
  for (unsigned int k=0;k<segments.size();k++){
    for (unsigned int n=0;n<segments[k]->hits.size();n++){
      const DFDCPseudo *fdchit=segments[k]->hits[n];
      fit.AddHit(fdchit);
      //bfield->GetField(x,y,z,Bx,By,Bz);
      //Bz_avg-=Bz;
      Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),fdchit->wire->origin.z());
    }
    num_hits+=segments[k]->hits.size();
  }   
  Bz=fabs(Bz)/double(num_hits);
  
  // Fit the points to a circle
  if (fit.FitCircleRiemann(segments[0]->rc)==NOERROR){
    double p=0.003*Bz*fit.r0/cos(atan(fit.tanl));
  
    if (p>3.){ // momentum is suspiciously high for a track going through both
      // the FDC and the CDC... try alternate circle fit...  
      fit.FitCircle();
    }
    return NOERROR;
  }
  
  return RESOURCE_UNAVAILABLE;   
}

// Method to match FDC and CDC candidates that projects the CDC candidate to the
// cdc endplate, where the FDC candidates are reported.
bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan,
                  vector<unsigned int> &cdc_forward_ids,
                  vector<DVector3>&cdc_endplate_projections,
                  vector<int>&cdc_forward_matches
                  ){
  // Momentum and position vectors for the FDC candidate
  DVector3 mom=fdccan->momentum();
  DVector3 pos=fdccan->position();
  double p_fdc=mom.Mag();
  ProjectHelixToZ(cdc_endplate.z(),fdccan->charge(),mom,pos); 

  // loop over the cdc candidates looking for the smallest distance
  // between the cdc and fdc projections to the end plate
  double diff_min=1000.; // candidate matching difference
  unsigned int jmin=0;
  for (unsigned int j=0;j<cdc_forward_ids.size();j++){
    if (cdc_forward_matches[j]) continue;
    double diff=(cdc_endplate_projections[j]-pos).Mag();
    if (diff<diff_min){
      diff_min=diff;
      jmin=j;
    }
  }
  
  if (DEBUG_HISTS){
    match_dist->Fill(pos.Perp(),diff_min);
    match_dist_vs_p->Fill(p_fdc,diff_min);
  }

  // Magnitude of the momentum of the potential cdc match
  double p_cdc=cdctrackcandidates[cdc_forward_ids[jmin]]->momentum().Mag();
  if (cdc_fdc_match(p_fdc,p_cdc,diff_min)){
    // Get the segment data
    vector<const DFDCSegment *>segments;
    fdccan->GetT(segments);

    // The following code snippet is intended to prevent matching a cdc
    // track candidate to a low momentum particle curling from the cdc endplate
    if (segments.size()>1){
      double theta=180./M_PI*mom.Theta(); 
      if (p_fdc<0.3 && theta<5.) return false;
    }

    unsigned int cdc_index=cdc_forward_ids[jmin];
    if (MakeCandidateFromMethod1(mom.Theta(),segments,
             cdctrackcandidates[cdc_index])){
         
      if (DEBUG_LEVEL>0)
   _DBG_ << ".. matched to CDC candidate #" << cdc_index <<endl;   
      
      // mark the cdc candidate as matched
      cdc_forward_matches[jmin]=1;

      return true;
    }
  } // got a cdc-fdc match
  
  // no match
  return false;
}

// Create new candidate after using Match Method 1 to match cdc and fdc 
// candidates
bool DTrackCandidate_factory::MakeCandidateFromMethod1(double theta,vector<const DFDCSegment *>&segments,const DTrackCandidate *cdccan){
  // JANA does not maintain the order that the segments were added
  // as associated objects. Therefore, we need to reorder them here
  // so segment[0] is the most upstream segment.
  stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);
  
  // Get the associated cdc hits
  vector<const DCDCTrackHit *>cdchits;
  cdccan->GetT(cdchits);
  
  // Sort CDC hits by layer
  stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing);
  
  // Abort if the track would have to pass from one quadrant into another
  // quadrant (this is more likely two roughly back-to-back tracks rather than 
  // a single track).
  const DVector2 fdc_hit_pos=segments[0]->hits[0]->xy;
  const DVector3 cdc_wire_origin=cdchits[0]->wire->origin;  
  if (cdc_wire_origin.x()*fdc_hit_pos.X()<0.
      || cdc_wire_origin.y()*fdc_hit_pos.Y()<0.){
    if (DEBUG_LEVEL>0) _DBG_ << "Skipping match of potential back-to-back tracks." <<endl;
    return false; 
  }
  
  // average Bz
  double Bz_avg=0.;
  
  // Redo circle fit with additional hits
  DHelicalFit fit;
  if (DoRefit(fit,segments,cdchits,Bz_avg)==NOERROR){
    // Determine the polar angle
    double theta_cdc=cdccan->momentum().Theta();
    if (segments.size()==1 && theta_cdc<M_PI_4){
      double numcdc=double(cdchits.size());
      double numfdc=segments[0]->hits.size();
      theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
    }
    fit.tanl=tan(M_PI_2-theta);
    
    // FDC hit
    const DFDCPseudo *fdchit=segments[0]->hits[0];
    double zhit=fdchit->wire->origin.z();
    DVector3 pos(fdchit->xy.X(),fdchit->xy.Y(),zhit);
    DVector3 mom;
    UpdatePositionAndMomentum(fit,Bz_avg,cdchits[0]->wire->origin,pos,mom);

    // Create new track candidate object 
    DTrackCandidate *can = new DTrackCandidate;
    can->used_cdc_indexes=cdccan->used_cdc_indexes;
    // circle parameters
    can->rc=fit.r0;
    can->xc=fit.x0;
    can->yc=fit.y0;
    
    // Add cdc and fdc hits to the track as associated objects
    for (unsigned int m=0;m<segments.size();m++){
      for (unsigned int n=0;n<segments[m]->hits.size();n++){
   can->AddAssociatedObject(segments[m]->hits[n]);
      }
    }
    for (unsigned int n=0;n<cdchits.size();n++){
      can->AddAssociatedObject(cdchits[n]); 
    }

    can->chisq=fit.chisq;
    can->Ndof=fit.ndof;
    Particle_t locPID = (FactorForSenseOfRotation*fit.h > 0.0) ? PiPlus : PiMinus;
    can->setPID(locPID);
    can->setMomentum(mom);
    can->setPosition(pos);
    
    trackcandidates.push_back(can);       
    
    return true;
  }
  return false;
}
 
// Method to match CDC and FDC candidates that projects an FDC candidate into 
// the CDC region using a helical approximation to the trajectory.  The code
// computes a doca to each wire in the cdc candidate and counts matches based
// on the probability that the cdc wire is on the track.
 bool DTrackCandidate_factory::MatchMethod2(const DTrackCandidate *fdccan,
                   const DTrackCandidate *cdccan
                   ){
   // Get the cdc hits associated with this cdc candidate
   vector<const DCDCTrackHit *>cdchits;
   cdccan->GetT(cdchits);   
   stable_sort(cdchits.begin(),cdchits.end(),CDCHitSortByLayerincreasing);
   
   // Variables to count the number of matching hits and the match fraction
   unsigned int num_match=0;
   unsigned int num_cdc=0;
   
     // Get the track parameters from the fdc candidate
   DVector3 pos=fdccan->position();
   DVector3 mom=fdccan->momentum();
   double q=fdccan->charge();

   // Check to see if the outer hit in the CDC does not exceed the radius of  
   // the FDC position to try to avoid false matches...
   unsigned int outer_index=cdchits.size()-1;
   DVector3 origin=cdchits[outer_index]->wire->origin;
   DVector3 dir=(1./cdchits[outer_index]->wire->udir.z())*cdchits[outer_index]->wire->udir;
   DVector3 cdc_outer_wire_pos=origin+(167.-origin.z())*dir;
   if (cdc_outer_wire_pos.Perp()>pos.Perp()) {
     return false;
   }
   // Make sure that the rough position of the CDC track at the end plate 
   // is not too far off from the position of the fdc track near the end plate
   if ((cdc_outer_wire_pos-pos).Mag()>5.){
     return false;
   }

   // loop over the cdc hits and count hits that agree with a projection of 
   // the helix into the cdc 
   for (unsigned int m=0;m<cdchits.size();m++){
     double variance=1.6;
     // Use a helical approximation to the track to match both axial and
     // stereo wires
     double dr2=DocaToHelix(cdchits[m],q,pos,mom);
     double prob=isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;
     
     if (prob>0.01) num_match++;
     if (DEBUG_LEVEL>1) _DBG_ << "CDC s: " << cdchits[m]->wire->straw
               << " r: " << cdchits[m]->wire->ring
            << " prob: " << prob << endl;
     
     num_cdc++;
   }
   
   if (num_match>=3 && double(num_match)/double(num_cdc)>0.33){              
     // Get the segment data
     vector<const DFDCSegment *>segments;
     fdccan->GetT(segments);
     
     // JANA does not maintain the order that the segments were added
     // as associated objects. Therefore, we need to reorder them here
     // so segment[0] is the most upstream segment.
     stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);
     
     // Abort if the track would have to pass from one quadrant into the opposite 
     // quadrant (this is more likely two roughly back-to-back tracks rather than 
     // a single track).
     const DVector2 fdc_hit_pos=segments[0]->hits[0]->xy;
     const DVector3 cdc_wire_origin=cdchits[0]->wire->origin;  
     if (cdc_wire_origin.x()*fdc_hit_pos.X()<0.
    && cdc_wire_origin.y()*fdc_hit_pos.Y()<0.){
       if (DEBUG_LEVEL>0) _DBG_ << "Skipping match of potential back-to-back tracks." <<endl;
       return false;
     }
     
     // Put the fdc candidate in the combined list
     DTrackCandidate *can = new DTrackCandidate;
     can->setPID((q > 0.0) ? PiPlus : PiMinus);
     
     // copy the list of cdc indices over from the cdc candidate
     can->used_cdc_indexes=cdccan->used_cdc_indexes;
     
     // Add the fdc hits to track candidate as associated objects 
     unsigned int num_fdc_hits=0;
     for (unsigned int m=0;m<segments.size();m++){
       for (unsigned int n=0;n<segments[m]->hits.size();n++){
    can->AddAssociatedObject(segments[m]->hits[n]);
    num_fdc_hits++;
       }
     }
     // Add the CDC hits to the track candidate
     for (unsigned int m=0;m<cdchits.size();m++){
       can->AddAssociatedObject(cdchits[m]);
     }
       
     // Average Bz
     double Bz_avg=0.;
     
     // Redo circle fit with additional hits
     DHelicalFit fit;
     //...first add the tangent to the dip angle to the fit class...
     double theta=fdccan->momentum().Theta();
     double theta_cdc=cdccan->momentum().Theta();
     if (segments.size()==1&& theta_cdc<M_PI_4){
       double numcdc=double(cdchits.size());
       double numfdc=segments[0]->hits.size();
       theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
     }
     fit.tanl=tan(M_PI_2-theta);
     if (DoRefit(fit,segments,cdchits,Bz_avg)==NOERROR){
       fit.FindSenseOfRotation();
       Particle_t locPID = ((FactorForSenseOfRotation*fit.h > 0.0) ? PiPlus : PiMinus);
       can->setPID(locPID);

       // FDC hit
       const DFDCPseudo *fdchit=segments[0]->hits[0];
       double zhit=fdchit->wire->origin.z();
       pos.SetXYZ(fdchit->xy.X(),fdchit->xy.Y(),zhit);
       UpdatePositionAndMomentum(fit,Bz_avg,cdchits[0]->wire->origin,pos,mom);
       
       // circle parameters
       can->rc=fit.r0;
       can->xc=fit.x0;
       can->yc=fit.y0;
       
       can->chisq=fit.chisq;
       can->Ndof=fit.ndof;
       can->setMomentum(mom);
       can->setPosition(pos);
     }
     else{
       //_DBG_ << endl;
       // circle parameters
       can->rc=fdccan->rc;
       can->xc=fdccan->xc;
       can->yc=fdccan->yc;
       can->Ndof=fdccan->Ndof;
       can->chisq=fdccan->chisq;
       can->setMomentum(fdccan->momentum());
       can->setPosition(fdccan->position());     
     }
     
     trackcandidates.push_back(can);
 
     return true;
   }     

   return false;
 }
 
 // Method to match CDC and FDC candidates that projects the CDC track into the
 // FDC region.
 bool DTrackCandidate_factory::MatchMethod3(const DTrackCandidate *cdccan,
                   vector<int> &forward_matches
                   ){
   // Get hits already linked to this candidate from associated objects
   vector<const DCDCTrackHit *>cdchits;
   cdccan->GetT(cdchits);
   stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing);
   
   // loop over fdc candidates
   for (unsigned int k=0;k<fdctrackcandidates.size();k++){
     if (forward_matches[k]==0){
       const DTrackCandidate *fdccan = fdctrackcandidates[k];
       // fdc and cdc candidate momenta
       double p_fdc=fdccan->momentum().Mag();
       double p_cdc=cdccan->momentum().Mag();
       if (p_fdc/p_cdc>0.5){       
    // Get the segment data
    vector<const DFDCSegment *>segments;
    fdccan->GetT(segments);
    stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);

    // Only try to match candidate if this fdc candidate has hits in the 
    // first package
    if (segments[0]->package>0) continue;
    
    
    // Abort if the track would have to pass from one quadrant into the opposite 
    // quadrant (this is more likely two roughly back-to-back tracks rather than 
    // a single track).
    const DVector2 fdc_hit_pos=segments[0]->hits[0]->xy;
    const DVector3 cdc_wire_origin=cdchits[0]->wire->origin;  
    if (cdc_wire_origin.x()*fdc_hit_pos.X()<0.
        && cdc_wire_origin.y()*fdc_hit_pos.Y()<0.){
      if (DEBUG_LEVEL>0) _DBG_ << "Skipping match of potential back-to-back tracks." <<endl;
      continue; 
    }
       
    // Momentum and position vectors for the CDC candidate
    DVector3 mom=cdccan->momentum();
    DVector3 pos=cdccan->position();
    double q=cdccan->charge();

    // Try to match unmatched fdc candidates
    int num_match=0;
    int num_hits=0;
    for (unsigned int m=0;m<segments.size();m++){
      for (unsigned int n=0;n<segments[m]->hits.size();n++){
        unsigned int ind=segments[m]->hits.size()-1-n;
        const DFDCPseudo *hit=segments[m]->hits[ind];
        
        if (pos.Perp()<48.5){
          ProjectHelixToZ(hit->wire->origin.z(),q,mom,pos);
          
          // difference
          DVector2 XY=hit->xy;
          double dx=XY.X()-pos.x();
          double dy=XY.Y()-pos.y();
          double dr2=dx*dx+dy*dy;
          
          double variance=1.0;
          double prob = isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;
          if (prob>0.01) num_match++;
          
          num_hits++;
        }
      }
      if (num_match>=3 
          || (num_match>0 && double(num_match)/double(num_hits)>0.33)){
        DTrackCandidate *can = new DTrackCandidate;
        
        can->setPID(cdccan->PID());
        can->used_cdc_indexes=cdccan->used_cdc_indexes;
        
        // mark the fdc track candidate as matched
        forward_matches[k]=1; 
        
        // Add cdc hits to the candidate
        for (unsigned int m=0;m<cdchits.size();m++){
          can->AddAssociatedObject(cdchits[m]);
        } 
        
        // Add fdc hits to the candidate
        for (unsigned int m=0;m<segments.size();m++){
          for (unsigned int n=0;n<segments[m]->hits.size();n++){
       can->AddAssociatedObject(segments[m]->hits[n]);
          }
        }
        
        // average Bz
        double Bz_avg=0.;
        
        // Redo helical fit with all available hits
        DHelicalFit fit; 
        if (DoRefit(fit,segments,cdchits,Bz_avg)==NOERROR){
          // Determine the polar angle
          double theta=fdccan->momentum().Theta();
          double theta_cdc=cdccan->momentum().Theta();
          if (segments.size()==1 && theta_cdc<M_PI_4){
       double numcdc=double(cdchits.size());
       double numfdc=segments[0]->hits.size();
       theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
          }
          fit.tanl=tan(M_PI_2-theta);
          
          // Redo the line fit
          fit.FitLineRiemann();
          
          // Set the charge
          fit.FindSenseOfRotation();   
          Particle_t locPID = ((FactorForSenseOfRotation*fit.h > 0.0) ? PiPlus : PiMinus);
          can->setPID(locPID);
          // FDC hit
          const DFDCPseudo *fdchit=segments[0]->hits[0];
          double zhit=fdchit->wire->origin.z();
          pos.SetXYZ(fdchit->xy.X(),fdchit->xy.Y(),zhit); 
          UpdatePositionAndMomentum(fit,Bz_avg,cdchits[0]->wire->origin,
                pos,mom);

          // circle parameters
          can->rc=fit.r0;
          can->xc=fit.x0;
          can->yc=fit.y0;

          can->chisq=fit.chisq;
          can->Ndof=fit.ndof;
          can->setMomentum(mom);
          can->setPosition(pos);
        }
        else{
          //_DBG_ << endl;
            // circle parameters
          can->rc=cdccan->rc;
          can->xc=cdccan->xc;
          can->yc=cdccan->yc;
          can->Ndof=cdccan->Ndof;
          can->chisq=cdccan->chisq;
          can->setMomentum(cdccan->momentum());
          can->setPosition(cdccan->position());
        }
      
        trackcandidates.push_back(can);
        
        if (DEBUG_LEVEL>0) _DBG_ << ".. matched to FDC candidate #" << k <<endl;
      
        return true;
      } // got a match
    } // loop over segments
       } // check that we have not already matched to this fdc candidate
     } // momentum check
   } // loop over fdc candidates
   
   return false;
 }

// This method tries to match unmatched fdc candidates 
bool DTrackCandidate_factory::MatchMethod4(const DTrackCandidate *srccan, 
                  vector<int> &forward_matches,
                  int &num_fdc_cands_remaining){
  // Get segments already linked to this candidate from associated objects
  vector<const DFDCSegment *>src_segments;
  srccan->GetT(src_segments);

  if (src_segments.size()==1) return false;
  stable_sort(src_segments.begin(), src_segments.end(), SegmentSortByLayerincreasing);
  
  if (DEBUG_LEVEL>0) _DBG_ << "Attempting matching method #4..." <<endl; 

  unsigned int pack1_last=src_segments[src_segments.size()-1]->package;
  unsigned int pack1_first=src_segments[0]->package;
  
  for (unsigned int k=0;k<fdctrackcandidates.size();k++){
    if (num_fdc_cands_remaining==0) break;
    if (forward_matches[k]==0){
      const DTrackCandidate *fdccan = fdctrackcandidates[k];

      // Get the segment data
      vector<const DFDCSegment *>segments;
      fdccan->GetT(segments);
      stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);
   
      const DFDCPseudo *firsthit=segments[0]->hits[0];
      unsigned int pack2_first=segments[0]->package;
      unsigned int pack2_last=segments[segments.size()-1]->package;
      
      // if (pack2_first>pack1_last || pack2_last<pack1_first){
      if (pack2_first-pack1_last==1 || pack1_first-pack2_last==1){
   // Momentum and position vectors for the input candidate
   DVector3 mom=srccan->momentum();
   DVector3 pos=srccan->position();
   double q=srccan->charge();
   DVector3 norm(0.,0.,1.);

   if (pack2_last<pack1_first){
     mom=(-1.0)*mom;
   }

   // Swim to the first hit in the candidate to which we wish to match
   // using the stepper
   stepper->SetCharge(q);
   stepper->SwimToPlane(pos,mom,firsthit->wire->origin,norm,NULL);
     
   double dx=pos.x()-firsthit->xy.X();
   double dy=pos.y()-firsthit->xy.Y();
   double d2=dx*dx+dy*dy;
   double variance=1.;
   double prob=TMath::Prob(d2/variance,1);
     
   unsigned int num_match=(prob>0.01)?1:0;

   // Try to match unmatched fdc candidates
   for (unsigned int i=1;i<segments[0]->hits.size();i++){
     const DFDCPseudo *hit=segments[0]->hits[i];
     
     if (pos.Perp()<48.5){
       ProjectHelixToZ(hit->wire->origin.z(),q,mom,pos);
       
       DVector2 XY=hit->xy;
       double dx=XY.X()-pos.x();
       double dy=XY.Y()-pos.y();
       double dr2=dx*dx+dy*dy;
       
       double variance=1.0;
       double prob = isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;
       if (prob>0.01) num_match++;
     }
   } 
   if (num_match>=3){
     forward_matches[k]=1;
       
     // Create new candidate for output
     DTrackCandidate *can = new DTrackCandidate;
   
     // average Bz
     double Bz_avg=0.;
     unsigned int num_hits=0;
       
     // Create a new DHelicalFit object for fitting combined data
     DHelicalFit fit; 
     // Fake point at origin
     fit.AddHitXYZ(0.,0.,TARGET_Z,BEAM_VAR,BEAM_VAR,0.,true); 
     // Add hits to the fit object and also to the track candidate
     // itself as associated objects
     for (unsigned int m=0;m<src_segments.size();m++){
       for (unsigned int n=0;n<src_segments[m]->hits.size();n++){
         const DFDCPseudo *fdchit=src_segments[m]->hits[n];
         fit.AddHit(fdchit);
         can->AddAssociatedObject(fdchit);
         
         //bfield->GetField(x,y,z,Bx,By,Bz);
         //Bz_avg-=Bz;
         Bz_avg+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),
                fdchit->wire->origin.z());
       }
       num_hits+=src_segments[m]->hits.size();
     }
     for (unsigned int m=0;m<segments.size();m++){
       for (unsigned int n=0;n<segments[m]->hits.size();n++){
         const DFDCPseudo *fdchit=segments[m]->hits[n];
         fit.AddHit(fdchit);
         can->AddAssociatedObject(fdchit);
         
         //bfield->GetField(x,y,z,Bx,By,Bz);
         //Bz_avg-=Bz;
         Bz_avg+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),
                      fdchit->wire->origin.z());
       }
       num_hits+=segments[m]->hits.size();
     }
      
     // Fit the points to a circle
     if (fit.FitCircleRiemann(segments[0]->rc)==NOERROR){
       // Compute new transverse momentum
       Bz_avg=fabs(Bz_avg)/double(num_hits);
         
       // Guess for theta and z from input candidates
       double theta=fdccan->momentum().Theta();  
       fit.tanl=tan(M_PI_2-theta);
       fit.z_vertex=srccan->position().Z();
       
       // Redo line fit
       fit.FitLineRiemann();

       double p=0.003*fit.r0*Bz_avg/cos(atan(fit.tanl));
       if (p>10.){ // Check for extremely stiff tracks, some of which 
         // will have unphysical momenta... try alternate circle fit 
         fit.FitCircle();
       }
       // Guess charge from fit
       fit.h=GetSenseOfRotation(fit,firsthit,srccan->position());
       Particle_t locPID = ((FactorForSenseOfRotation*fit.h > 0.0) ? PiPlus : PiMinus);
      can->setPID(locPID);

       // put z position just upstream of the first hit in z
       const DHFHit_t *myhit=fit.GetHits()[0];
       pos.SetXYZ(myhit->x,myhit->y,myhit->z);
       GetPositionAndMomentum(myhit->z-1.,fit,Bz_avg,pos,mom);

       // circle parameters
       can->rc=fit.r0;
       can->xc=fit.x0;
       can->yc=fit.y0;
       
       can->setPosition(pos);
       can->setMomentum(mom);
       can->chisq=fit.chisq;
       can->Ndof=fit.ndof;
       trackcandidates.push_back(can);
       
       num_fdc_cands_remaining--;
       
       if (DEBUG_LEVEL>0)
         _DBG_ << "Found a match using method #4" <<endl;
       return true;
     } // circle fit
   } // got a match?
      } // is the first package in the new fdc candidate downstream of the last?
    } // check that we have not already matched this fdc candidate
  } // loop over fdc track candidates

  return false;
}

// Matching method to try to deal with CDC candidates that do not point 
// directly at the the cdc endplate but could possibly be matched to FDC 
// candidates.
bool DTrackCandidate_factory::MatchMethod5(DTrackCandidate *can,
                  vector<const DCDCTrackHit *>&cdchits,
                  vector<int> &forward_matches
                  ){
  // Loop over the fdc track candidates
  for (unsigned int k=0;k<fdctrackcandidates.size();k++){
    if (forward_matches[k]==0){
      const DTrackCandidate *fdccan = fdctrackcandidates[k];
      
      unsigned int num_match=0;
      unsigned int num_cdc=0;
      
      DVector3 pos=fdccan->position();
      DVector3 mom=fdccan->momentum();
      double q=fdccan->charge();

      for (unsigned int m=0;m<cdchits.size();m++){
   double variance=1.0;
   // Use a helical approximation to the track to match both axial and
   // stereo wires
   double dr2=DocaToHelix(cdchits[m],q,pos,mom);
   double prob=isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;
   
   if (prob>0.01) num_match++;
   if (DEBUG_LEVEL>1) _DBG_ << "CDC s: " << cdchits[m]->wire->straw
             << " r: " << cdchits[m]->wire->ring
             << " prob: " << prob << endl;

   num_cdc++;
      }

      if (num_match>=3 && double(num_match)/double(num_cdc)>0.33){   
   // Get the segment data
   vector<const DFDCSegment *>segments;
   fdccan->GetT(segments);
   
   // JANA does not maintain the order that the segments were added
   // as associated objects. Therefore, we need to reorder them here
   // so segment[0] is the most upstream segment.
   stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);

   // Abort if the track would have to pass from one quadrant into the opposite 
   // quadrant (this is more likely two roughly back-to-back tracks rather than 
   // a single track).
   const DVector2 fdc_hit_pos=segments[0]->hits[0]->xy;
   const DVector3 cdc_wire_origin=cdchits[0]->wire->origin;  
   if (cdc_wire_origin.x()*fdc_hit_pos.X()<0.
       && cdc_wire_origin.y()*fdc_hit_pos.Y()<0.){
     if (DEBUG_LEVEL>0) _DBG_ << "Skipping match of potential back-to-back tracks." <<endl;
     continue; 
   }

   // Mark this fdc candidate as having a match to a cdc candidate
   forward_matches[k]=1;
   
   if (DEBUG_LEVEL>0) _DBG_ << "... matched to FDC candidate #" << k <<endl;
   
   // Add fdc hits to the candidate
   for (unsigned int m=0;m<segments.size();m++){
     for (unsigned int n=0;n<segments[m]->hits.size();n++){
       can->AddAssociatedObject(segments[m]->hits[n]);
     }
   }
   
   // Average Bz
   double Bz_avg=0.;
   
   // Instantiate the helical fitter to do the refit
   DHelicalFit fit;
   if (DoRefit(fit,segments,cdchits,Bz_avg)==NOERROR){
     // circle parameters
     can->rc=fit.r0;
     can->xc=fit.x0;
     can->yc=fit.y0;
     
     // Determine the polar angle
     double theta=fdccan->momentum().Theta();
     fit.tanl=tan(M_PI_2-theta);
     
     Particle_t locPID = ((FactorForSenseOfRotation*fit.h > 0.0) ? PiPlus : PiMinus);
     can->setPID(locPID);

     // FDC hit
     const DFDCPseudo *fdchit=segments[0]->hits[0];
     double zhit=fdchit->wire->origin.z();
     pos.SetXYZ(fdchit->xy.X(),fdchit->xy.Y(),zhit);
     UpdatePositionAndMomentum(fit,Bz_avg,cdchits[0]->wire->origin,
                pos,mom);

     can->chisq=fit.chisq;
     can->Ndof=fit.ndof;
     can->setMomentum(mom);
     can->setPosition(pos);
   }
   return true;
      } // check for match
    } // check that the candidate has not already been matched
  } // loop over fdc candidates
  
  return false;
}

// Method to match FDC candidates with stray CDC hits unassociated with track
// candidates
void DTrackCandidate_factory::MatchMethod6(DTrackCandidate *can, 
                  vector<const DFDCPseudo*>&fdchits,
                  vector<unsigned int>&used_cdc_hits,  
                  unsigned int &num_unmatched_cdcs
                  ){
  // Input candidate parameters
  DVector3 pos=can->position();
  DVector3 mom=can->momentum();
  double q=can->charge();

 // Variables to keep track of cdc hits 
  bool got_inner_index=false;
  unsigned int inner_index=0;
  int id_for_smallest_dr=-1;
  int old_ring=-1;
  double dr2_min=1e6;
  
  for (unsigned int k=0;k<used_cdc_hits.size();k++){
    // We only want to use one hit from a given ring to avoid confusing the 
    // circle refit with curlers...
    if (mycdchits[k]->wire->ring!=old_ring){
      if (id_for_smallest_dr>=0){
   if (!used_cdc_hits[id_for_smallest_dr]){
     num_unmatched_cdcs--;
     used_cdc_hits[id_for_smallest_dr]=1;
     
     can->used_cdc_indexes.push_back(id_for_smallest_dr);
     can->AddAssociatedObject(mycdchits[id_for_smallest_dr]);
     
     if (got_inner_index==false){
       inner_index=id_for_smallest_dr;
       got_inner_index=true;
     }
   }
      }
      // Reset matching variables
      dr2_min=1e6;
      id_for_smallest_dr=-1;
    }
    
    if (!used_cdc_hits[k]){
      double variance=1.0;
      // Use a helical approximation to the track to match both axial and
      // stereo wires
      double dr2=DocaToHelix(mycdchits[k],q,pos,mom);
      double prob=isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;

      if (prob>0.5){
   // Look for closest match
   if (dr2<dr2_min){
     dr2_min=dr2;
     id_for_smallest_dr=k;
   }
      }
    }  
    old_ring=mycdchits[k]->wire->ring;
  } 
  // Grab the last potential match
  if (id_for_smallest_dr>=0){
    if (!used_cdc_hits[id_for_smallest_dr]){
      num_unmatched_cdcs--;
      used_cdc_hits[id_for_smallest_dr]=1;
      
      can->used_cdc_indexes.push_back(id_for_smallest_dr);
      can->AddAssociatedObject(mycdchits[id_for_smallest_dr]);
    }
  }

  // We found some matched hits.  Update the start position of the candidate.
  if (got_inner_index){
    const DFDCPseudo *firsthit=fdchits[0];

    // Compute the average Bz
    double Bz_avg=0.,denom=0.;
    for (unsigned int m=0;m<fdchits.size();m++){
      const DFDCPseudo *fdchit=fdchits[m];
      double x=fdchit->xy.X();
      double y=fdchit->xy.Y();
      double z=fdchit->wire->origin.z();
      
      Bz_avg+=bfield->GetBz(x,y,z);
      denom+=1.;
    }
    Bz_avg=fabs(Bz_avg)/denom;

    // Store the current track parameters in the DHelicalFit class
    DHelicalFit fit;
    fit.r0=mom.Perp()/(0.003*Bz_avg);
    fit.phi=mom.Phi();
    fit.h=q*FactorForSenseOfRotation;
    fit.x0=pos.x()-fit.h*fit.r0*sin(fit.phi);
    fit.y0=pos.y()+fit.h*fit.r0*cos(fit.phi);
    fit.tanl=tan(M_PI_2-mom.Theta());
    fit.z_vertex=0; // actualy we don't need this..
    // Use fdc hit for input to the following routines because the z-position is 
    // well-defined...
    double zhit=firsthit->wire->origin.z();
    pos.SetXYZ(firsthit->xy.X(),firsthit->xy.Y(),zhit);
    UpdatePositionAndMomentum(fit,Bz_avg,mycdchits[inner_index]->wire->origin,
               pos,mom);

    // update the track parameters
    can->setMomentum(mom);
    can->setPosition(pos);
    
    if (DEBUG_LEVEL>0) _DBG_ << "... matched stray CDC hits ..." << endl;
  } // match at least one cdc hit
}

// This method tries to match unmatched fdc candidates to candidates that
// already have fdc/cdc matches 
bool DTrackCandidate_factory::MatchMethod7(DTrackCandidate *srccan, 
                  vector<int> &forward_matches,
                  int &num_fdc_cands_remaining){ 
  if (DEBUG_LEVEL>0)  _DBG_ << "Attempting matching method #7..." <<endl; 

  // Get the hits associated with this candidate
  vector<const DFDCPseudo *>fdchits;
  srccan->GetT(fdchits);
  if (fdchits.size()==0) return false;

  stable_sort(fdchits.begin(),fdchits.end(),FDCHitSortByLayerincreasing);
  unsigned int pack1_last=(fdchits[fdchits.size()-1]->wire->layer-1)/6;
  
  for (unsigned int k=0;k<fdctrackcandidates.size();k++){
    if (num_fdc_cands_remaining==0) break;
    if (forward_matches[k]==0){
      const DTrackCandidate *fdccan = fdctrackcandidates[k];

      // Get the segment data
      vector<const DFDCSegment *>segments;
      fdccan->GetT(segments);
      stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);
   
      const DFDCPseudo *firsthit=segments[0]->hits[0];
      unsigned int pack2_first=segments[0]->package;
      
      if (pack2_first-pack1_last==1){ // match in adjacent packages
   // Momentum and position vectors for the input candidate
   DVector3 mom=srccan->momentum();
   DVector3 pos=srccan->position();
   double q=srccan->charge();
   DVector3 norm(0.,0.,1.);

   // Swim to the first hit in the candidate to which we wish to match
   // using the stepper
   stepper->SetCharge(q);
   stepper->SwimToPlane(pos,mom,firsthit->wire->origin,norm,NULL);
     
   double dx=pos.x()-firsthit->xy.X();
   double dy=pos.y()-firsthit->xy.Y();
   double d2=dx*dx+dy*dy;
   double variance=1.;
   double prob=TMath::Prob(d2/variance,1);
     
   unsigned int num_match=(prob>0.01)?1:0;

   // Try to match unmatched fdc candidates
   for (unsigned int i=1;i<segments[0]->hits.size();i++){
     const DFDCPseudo *hit=segments[0]->hits[i];
     
     if (pos.Perp()<48.5){
       ProjectHelixToZ(hit->wire->origin.z(),q,mom,pos);
       
       DVector2 XY=hit->xy;
       double dx=XY.X()-pos.x();
       double dy=XY.Y()-pos.y();
       double dr2=dx*dx+dy*dy;
       
       double variance=1.0;
       double prob = isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;
       if (prob>0.01) num_match++;
     }
   } 
   if (num_match>=3){
     forward_matches[k]=1;
     num_fdc_cands_remaining--;

     // average Bz
     double Bz=0.;
     unsigned int num_hits=0;
       
     // Create a new DHelicalFit object for fitting combined data
     DHelicalFit fit; 

     // Get the cdc hits associated with this track candidate and add the
     // axial straws to the list of hits to use in the circle fit
     vector<const DCDCTrackHit *>cdchits;
     srccan->GetT(cdchits);
     stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing);
     for (unsigned int i=0;i<cdchits.size();i++){
       if (cdchits[i]->is_stereo==false){
         double cov=0.213;  //guess
         const DVector3 origin=cdchits[i]->wire->origin;
         double x=origin.x(),y=origin.y(),z=origin.z();
         fit.AddHitXYZ(x,y,z,cov,cov,0.,true);
         Bz+=bfield->GetBz(x,y,z);
         num_hits++;
       }   
     }     
     // Add the FDC hits from the existing track to this list
     for (unsigned int i=0;i<fdchits.size();i++){
       fit.AddHit(fdchits[i]);
       double zhit=fdchits[i]->wire->origin.z();
      
       Bz+=bfield->GetBz(fdchits[i]->xy.X(),fdchits[i]->xy.Y(),zhit);
       num_hits++;

       if (zhit<firsthit->wire->origin.z()){
         firsthit=fdchits[i];
       }
     }
     // Add the new set of FDC hits; also add them as associated objects
     // to the track
     for (unsigned int m=0;m<segments.size();m++){
       for (unsigned int n=0;n<segments[m]->hits.size();n++){
         const DFDCPseudo *fdchit=segments[m]->hits[n];
         fit.AddHit(fdchit);
         double zhit=fdchit->wire->origin.z();

         Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),zhit);
         srccan->AddAssociatedObject(fdchit);
       }
       num_hits+=segments[m]->hits.size();
     }   
     Bz=fabs(Bz)/double(num_hits);

     // Redo the circle fit with the extra hits
     if (fit.FitCircleRiemann(srccan->rc)==NOERROR){
       // Determine the polar angle
       double theta=srccan->momentum().Theta();
       fit.tanl=tan(M_PI_2-theta);

       double p=0.003*fit.r0*Bz/cos(atan(fit.tanl));
       if ((cdchits.size())>0?(p>3.):(p>10.)){
         // Suspiciously high momentum:  try alternate circle fit...
         fit.FitCircle();
       }    
      
       if (cdchits.size()>0){
         // FDC hit
         double zhit=firsthit->wire->origin.z();
         pos.SetXYZ(firsthit->xy.X(),firsthit->xy.Y(),zhit); 
         UpdatePositionAndMomentum(fit,Bz,cdchits[0]->wire->origin,
               pos,mom);
       }
       else{
         // put z position just upstream of the first hit in z
         const DHFHit_t *myhit=fit.GetHits()[0];
         pos.SetXYZ(myhit->x,myhit->y,myhit->z);
         GetPositionAndMomentum(myhit->z-1.,fit,Bz,pos,mom);
       }
       srccan->rc=fit.r0;
       srccan->xc=fit.x0;
       srccan->yc=fit.y0;
       
       srccan->chisq=fit.chisq;
       srccan->Ndof=fit.ndof;
       srccan->setMomentum(mom);
       srccan->setPosition(pos);
     }
       
     if (DEBUG_LEVEL>0) _DBG_ << "Found a match using method #7" <<endl;
     return true;
   } // got a match?
      } // is the first package in the new fdc candidate downstream of the last?
    } // check that we have not already matched this fdc candidate
  } // loop over fdc track candidates

  return false;
}


// This method tries to improve the quality of the circle fit of a cdc
// candidate by assuming that if the outermost hit is in a stereo straw, 
// if this track is going to match to an fdc candidate, the location of the 
// hit along the straw is likely to be near the end of the straw, thus 
// providing another useful point on the circle.  After refitting the circle,
// the code adjusts the dip angle and tries to match to the remaining fdc 
// candidates.
bool DTrackCandidate_factory::MatchMethod8(const DTrackCandidate *cdccan,
                  vector<int> &forward_matches
                  ){
  // Get the cdc hits associated with this track candidate    
  vector<const DCDCTrackHit *>cdchits;
  cdccan->GetT(cdchits);
  stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing);
  
  unsigned int last_index=cdchits.size()-1;
  if (cdchits[last_index]->is_stereo==true 
      && cdchits[last_index]->wire->ring<13){
    DHelicalFit fit;
    
    // Assume the outermost stereo hit occurs near the end of the straw
    DVector3 origin=cdchits[last_index]->wire->origin;
    DVector3 dir=cdchits[last_index]->wire->udir;
    DVector3 wirepos=origin+(75./dir.z())*dir;
    double cov=0.2;
    double x=wirepos.x(),y=wirepos.y(),z=wirepos.z();
    fit.AddHitXYZ(x,y,z,cov,cov,0,true);
    double Bz=bfield->GetBz(x,y,z);
    int num_hits=1;

    // Add the cdc axial wires to the list of hits to use in the fit 
    for (unsigned int k=0;k<cdchits.size();k++){   
      if (cdchits[k]->is_stereo==false){
   cov=0.2;  //guess
   origin=cdchits[k]->wire->origin;
   double x=origin.x(),y=origin.y(),z=origin.z();
   fit.AddHitXYZ(x,y,z,cov,cov,0.,true);
   Bz+=bfield->GetBz(x,y,z);
   num_hits++;
      }   
    }
    Bz=fabs(Bz)/num_hits;
    
    // Fit the points to a circle assuming that the circle goes through the 
    // origin.
    if (fit.FitCircle()!=NOERROR) return false;
     
    // Determine the tangent of the dip angle
    double tworc=2.*fit.r0;
    double ratio=wirepos.Perp()/tworc;
    double sperp=tworc*((ratio<1.)?asin(ratio):M_PI_2);
    fit.tanl=(wirepos.z()-TARGET_Z)/sperp;

    // Find sense of rotation (proportional to the charge)
    fit.GuessChargeFromCircleFit();
    double q=fit.h*FactorForSenseOfRotation;

    // Provide an initial guess for the position and momentum that is just 
    // outside the CDC tracking volume
    DVector3 pos=wirepos,mom;
    UpdatePositionAndMomentum(fit,Bz,cdchits[0]->wire->origin,pos,mom);

    // Set the charge for the stepper 
    stepper->SetCharge(q);
    // Vector normal to the fdc planes
    DVector3 norm(0.,0.,1.);
    
    // Loop over the fdc candidates looking for a match
    for (unsigned int k=0;k<fdctrackcandidates.size();k++){
      if (forward_matches[k]==0){
   const DTrackCandidate *fdccan = fdctrackcandidates[k];
     
   // Get the segment data
   vector<const DFDCSegment *>segments;
   fdccan->GetT(segments);
   stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);
   
   // only do this for candidates with segments in the first package
   if (segments[0]->package>0) continue;
   
   const DFDCPseudo *firsthit=segments[0]->hits[0];
   
   // Make copies of the momentum and position vectors for input to 
   // the swimmer
   DVector3 my_mom(mom);
   DVector3 my_pos(pos);
    
   // Swim to the first hit in the candidate to which we wish to match
   // using the stepper
   stepper->SwimToPlane(my_pos,my_mom,firsthit->wire->origin,norm,NULL);
   
   // Match based on proximity of projection to hit
   double dx=my_pos.x()-firsthit->xy.X();
   double dy=my_pos.y()-firsthit->xy.Y();
   double d2=dx*dx+dy*dy;
   double variance=1.;
   double prob=TMath::Prob(d2/variance,1);
   
   unsigned int num_match=(prob>0.01)?1:0;
   
   // Try to match further hits in most upstream segment
   for (unsigned int i=1;i<segments[0]->hits.size();i++){
     const DFDCPseudo *hit=segments[0]->hits[i];
     
     if (my_pos.Perp()<48.5){
       ProjectHelixToZ(hit->wire->origin.z(),q,my_mom,my_pos);
       
       DVector2 XY=hit->xy;
       double dx=XY.X()-my_pos.x();
       double dy=XY.Y()-my_pos.y();
       double dr2=dx*dx+dy*dy;
       
       double variance=1.0;
       double prob = isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;
       if (prob>0.01) num_match++;
     }
   } 
   if (num_match>=3){
     forward_matches[k]=1;
     
     // Keep track of magnetic field in FDC
     double Bz_fdc=0;
     unsigned int num_hits_fdc=0;
     
     // Drop the first cdc hit in the list, since it isn't an
     // axial straw so the actual x,y position of the wire 
     // depends on z, which we don't really know yet.
     fit.PruneHit(0);
     
     // Add the fdc hits to the fit class
     for (unsigned int m=0;m<segments.size();m++){
       for (unsigned int n=0;n<segments[m]->hits.size();n++){
         const DFDCPseudo *fdchit=segments[m]->hits[n];
         fit.AddHit(fdchit);
         
         Bz_fdc+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),
                fdchit->wire->origin.z());
       }
       num_hits_fdc+=segments[m]->hits.size();
     }
     
     // Fake point at origin
     fit.AddHitXYZ(0.,0.,TARGET_Z,BEAM_VAR,BEAM_VAR,0.,true);
     // Fit the points to a circle
     if (fit.FitCircleRiemann(segments[0]->rc)==NOERROR){
       // Use the fdc track candidate to get tanl
       double theta=fdccan->momentum().Theta();
       double theta_cdc=cdccan->momentum().Theta();
       if (segments.size()==1 && theta_cdc<M_PI_4){
         double numcdc=double(cdchits.size());
         double numfdc=segments[0]->hits.size();
         theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
       }
       fit.tanl=tan(M_PI_2-theta);
       
       Bz=0.5*(Bz+fabs(Bz_fdc)/num_hits_fdc);
       
       double p=0.003*fit.r0*Bz/cos(atan(fit.tanl));
       if (p>10.){
         // momentum is suspiciously high for a track going through both
         // the FDC and the CDC... try alternate circle fit
         fit.FitCircle();
       }
       // FDC hit
       const DFDCPseudo *fdchit=segments[0]->hits[0];
       double zhit=fdchit->wire->origin.z();
       pos.SetXYZ(fdchit->xy.X(),fdchit->xy.Y(),zhit);
       UpdatePositionAndMomentum(fit,Bz,cdchits[0]->wire->origin,pos,mom);
       
       DTrackCandidate *can = new DTrackCandidate;
       // circle parameters
       can->rc=fit.r0;
       can->xc=fit.x0;
       can->yc=fit.y0;
       
       can->setPID((q > 0.0) ? PiPlus : PiMinus);
       can->setMomentum(my_mom);
       can->setPosition(my_pos);
       can->chisq=fit.chisq;
       can->Ndof=fit.ndof;
       
       // Add hits as associated objects
       for (unsigned int m=0;m<segments.size();m++){
         for (unsigned int n=0;n<segments[m]->hits.size();n++){
      const DFDCPseudo *fdchit=segments[m]->hits[n];
      can->AddAssociatedObject(fdchit);
         }
       }
       for (unsigned int n=0;n<cdchits.size();n++){
         can->AddAssociatedObject(cdchits[n]);
       }
       
       trackcandidates.push_back(can);         
       
       if (DEBUG_LEVEL>0)  _DBG_ << "Matched using Method #8" <<endl;
       
       return true;
     } // circle fit
   } // match criterion
      } // check that fdc candidate has not already been matched
    } // loop over fdc candidates
  } // Check that the outermost hit is a stereo hit

  return false;
}


// Method to try to match unmatched FDC segments by forcing the circle fits 
// to go through the origin.  
bool DTrackCandidate_factory::MatchMethod9(unsigned int src_index,
                  const DTrackCandidate *srccan,
                  const DFDCSegment *segment,
                  vector<const DTrackCandidate*>&cands,
                  vector<int> &forward_matches){
  if (DEBUG_LEVEL>0){
    _DBG_ << "Attempting Match method #9..." << endl;
  }
  double q=srccan->charge();
  int pack1=segment->package;

  // Get hits from segment and redo fit forcing the circle to go through (0,0)
  DHelicalFit fit1;
  for (unsigned int n=0;n<segment->hits.size();n++){
    const DFDCPseudo *hit=segment->hits[n];
    fit1.AddHit(hit);
  }
  if (fit1.FitCircle()!=NOERROR) return false;

   // Sense of rotation
  fit1.h=q*FactorForSenseOfRotation;
  
  // Use the fdc track candidate to get tanl
  double theta=srccan->momentum().Theta();
  fit1.tanl=tan(M_PI_2-theta);

  // Find the magnetic field at the first hit in the first segment
  double x=segment->hits[0]->xy.X();
  double y=segment->hits[0]->xy.Y();
  double z=segment->hits[0]->wire->origin.z();
  double Bz=fabs(bfield->GetBz(x,y,z));

  // Get position and momentum for the first segment
  DVector3 mypos,mymom;
  GetPositionAndMomentum(z,fit1,Bz,mypos,mymom);

  // Loop over the rest of the fdc track candidates, skipping those that have
  // already been used
  for (unsigned int k=src_index+1;k<forward_matches.size();k++){
    if (forward_matches[k]==0){
      const DTrackCandidate *can2=cands[k];
      // Get the segment data
      vector<const DFDCSegment *>segments2;
      can2->GetT(segments2);
       
      int pack2=segments2[0]->package;
      if (abs(pack1-pack2)>0){
   // Get hits from the second segment and redo fit forcing circle 
   // to go through (0,0)
   DHelicalFit fit2;
   for (unsigned int n=0;n<segments2[0]->hits.size();n++){
     const DFDCPseudo *hit=segments2[0]->hits[n];
     fit2.AddHit(hit);
   }
   if (fit2.FitCircle()!=NOERROR) continue;
         
   // Match using centers of circles
   double dx=fit1.x0-fit2.x0;
   double dy=fit1.y0-fit2.y0;
   double circle_center_diff2=dx*dx+dy*dy;
   double got_match=false;
   if (circle_center_diff2<9.0) got_match=true;
   // try another matching method here if got_match==false 
   else{    
     // Sense of rotation
     double q2=can2->charge();
     fit2.h=q2*FactorForSenseOfRotation;
     
     // Use the fdc track candidate to get tanl
     theta=can2->momentum().Theta();
     fit2.tanl=tan(M_PI_2-theta);

     // Try to match segments by swimming through the field
     got_match=MatchMethod11(q,mypos,mymom,fit2,segment,segments2[0]);
     //    _DBG_ << endl;
   }
   if (got_match){
     forward_matches[k]=1;
     forward_matches[src_index]=1;
     
     // Create a new DTrackCandidate for output
     DTrackCandidate *can = new DTrackCandidate;
     
     // variables for finding <Bz>
     double Bz=0;
     int num_hits=0;
     
     // Add hits from first segment as associated objects
     for (unsigned int n=0;n<segment->hits.size();n++){
       const DFDCPseudo *fdchit=segment->hits[n];
       can->AddAssociatedObject(fdchit);
       
       Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),
               fdchit->wire->origin.z());
       num_hits++;
     }
     
     // Add the hits from the second segment to the first fit object
     // and refit the circle.  Also add them as associated objects
     // to the new candidate.
     for (unsigned int n=0;n<segments2[0]->hits.size();n++){
       const DFDCPseudo *hit=segments2[0]->hits[n];
       fit1.AddHit(hit); 
       can->AddAssociatedObject(hit);
   
       Bz+=bfield->GetBz(hit->xy.X(),hit->xy.Y(),
               hit->wire->origin.z());
       num_hits++;        
     }
     Bz=fabs(Bz)/double(num_hits);
     
     // Initialize variables needed for output
     DVector3 mom=srccan->momentum();
     DVector3 pos=srccan->position();

     can->Ndof=srccan->Ndof;
     can->chisq=srccan->chisq;
     
     if (fit1.FitCircleRiemann(fit1.r0)==NOERROR){    
       can->Ndof=fit1.ndof;
       can->chisq=fit1.chisq;
       
       // Redo line fit
       fit1.FitLineRiemann();
       
       double p=0.003*fit1.r0*Bz/cos(atan(fit1.tanl));
       if (p>10.){
         // Try an alternate circle fit
         fit1.FitCircle();
       }

       // Guess charge from fit
       fit1.h=GetSenseOfRotation(fit1,segments2[0]->hits[0],
                  srccan->position());
       q=FactorForSenseOfRotation*fit1.h;
       
       // put z position just upstream of the first hit in z
       const DHFHit_t *myhit=fit1.GetHits()[0];
       pos.SetXYZ(myhit->x,myhit->y,myhit->z);
       GetPositionAndMomentum(myhit->z-1.,fit1,Bz,pos,mom);
     }
     // circle parameters
     can->rc=fit1.r0;
     can->xc=fit1.x0;
     can->yc=fit1.y0;
     
     can->setPID((q > 0.0) ? PiPlus : PiMinus);
     can->setPosition(pos);
     can->setMomentum(mom);
     
     trackcandidates.push_back(can);
    
     if (DEBUG_LEVEL>0){
       _DBG_ << "Match method #9 succeeded..." << endl;
     }

     return true;
   } // circle center match
      } // different packages?
    } // already matched?
  } // loop over tracks

  return false;
}

// Method to try to match unmatched FDC segments by assuming that the tracks
// are sufficiently stiff we are better served by something closer to a
// "straight-line" approximation
bool DTrackCandidate_factory::MatchMethod10(unsigned int src_index,
                  const DTrackCandidate *srccan,
                  const DFDCSegment *segment,
                  vector<const DTrackCandidate*>&cands,
                  vector<int> &forward_matches){
  if (DEBUG_LEVEL>0){
    _DBG_ << "Attempting Match method #10..." << endl;
  }
  double q=srccan->charge();
  int pack1=segment->package;

  // Get hits from segment and redo fit
  DHelicalFit fit1;
  for (unsigned int n=0;n<segment->hits.size();n++){
    const DFDCPseudo *hit=segment->hits[n];
    fit1.AddHit(hit);
  }
  fit1.FitCircleStraightTrack();

  // Sense of rotation
  fit1.h=q*FactorForSenseOfRotation;
  
  // Use the fdc track candidate to get tanl
  double theta=srccan->momentum().Theta();
  fit1.tanl=tan(M_PI_2-theta);

  // Find the magnetic field at the first hit in the first segment
  double x=segment->hits[0]->xy.X();
  double y=segment->hits[0]->xy.Y();
  double z=segment->hits[0]->wire->origin.z();
  double Bz=fabs(bfield->GetBz(x,y,z));

  // Get position and momentum for the first segment
  DVector3 mypos,mymom;
  GetPositionAndMomentum(z,fit1,Bz,mypos,mymom);

  // Loop over the rest of the fdc track candidates, skipping those that have
  // already been used
  for (unsigned int k=src_index+1;k<forward_matches.size();k++){
    if (forward_matches[k]==0){
      const DTrackCandidate *can2=cands[k];
      // Get the segment data
      vector<const DFDCSegment *>segments2;
      can2->GetT(segments2);
       
      int pack2=segments2[0]->package;
      if (abs(pack1-pack2)>0){
   // Get hits from the second segment and redo fit
   DHelicalFit fit2;
   for (unsigned int n=0;n<segments2[0]->hits.size();n++){
     const DFDCPseudo *hit=segments2[0]->hits[n];
     fit2.AddHit(hit);
   }
   fit2.FitCircleStraightTrack();
   
   // Sense of rotation
   double q2=can2->charge();
   fit2.h=q2*FactorForSenseOfRotation;
   
   // Use the fdc track candidate to get tanl
   theta=can2->momentum().Theta();
   fit2.tanl=tan(M_PI_2-theta);

   // Try to match segments by swimming through the field
   if (MatchMethod11(q,mypos,mymom,fit2,segment,segments2[0])){
     forward_matches[k]=1;
     forward_matches[src_index]=1;

     // Create a new DTrackCandidate for output
     DTrackCandidate *can = new DTrackCandidate;
     
     // variables for finding <Bz>
     double Bz=0;
     int num_hits=0;
     
     // Add hits from first segment as associated objects
     for (unsigned int n=0;n<segment->hits.size();n++){
       const DFDCPseudo *fdchit=segment->hits[n];
       can->AddAssociatedObject(fdchit);
    
       Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),
               fdchit->wire->origin.z());
       num_hits++;
     }
     
     // Add the hits from the second segment to the first fit object
     // and refit the circle.  Also add them as associated objects
     // to the new candidate.
     for (unsigned int n=0;n<segments2[0]->hits.size();n++){
       const DFDCPseudo *hit=segments2[0]->hits[n];
       fit1.AddHit(hit); 
       can->AddAssociatedObject(hit);
     
       Bz+=bfield->GetBz(hit->xy.X(),hit->xy.Y(),
               hit->wire->origin.z());
       num_hits++;        
     }
     Bz=fabs(Bz)/double(num_hits);
     
     // Initialize variables needed for output
     DVector3 mom=srccan->momentum();
     DVector3 pos=srccan->position();
     double q=srccan->charge(); 
     can->Ndof=srccan->Ndof;
     can->chisq=srccan->chisq;
     
     if (fit1.FitCircleRiemann(fit1.r0)==NOERROR){
       can->Ndof=fit1.ndof;
       can->chisq=fit1.chisq;
         
       // Redo line fit
       fit1.FitLineRiemann();
         
       double p=0.003*fit1.r0*Bz/cos(atan(fit1.tanl));
       if (p>10.){
         // unphysical momentum, try alternate circle fit...
         fit1.FitCircle();
       }
       
       // Guess charge from fit
       fit1.h=GetSenseOfRotation(fit1,segments2[0]->hits[0],
                  srccan->position());
       q=FactorForSenseOfRotation*fit1.h;
       
       // put z position just upstream of the first hit in z
         const DHFHit_t *myhit=fit1.GetHits()[0];
         pos.SetXYZ(myhit->x,myhit->y,myhit->z);
         GetPositionAndMomentum(myhit->z-1.,fit1,Bz,pos,mom);
     }
     // circle parameters
     can->rc=fit1.r0;
     can->xc=fit1.x0;
     can->yc=fit1.y0;

     can->setPID((q > 0.0) ? PiPlus : PiMinus);
     can->setPosition(pos);
     can->setMomentum(mom);
     
     trackcandidates.push_back(can);

     if (DEBUG_LEVEL>0){
       _DBG_ << "Match method #10 succeeded..." << endl;
     }
     return true;
   } // minimum number of matching hits
      } // different packages?
    } // already matched?
  } // loop over tracks

  return false;
}

// Swims from one FDC segment to another segment looking for a match.  If this
// fails, swims from the second second to the first in the opposite direction. 
bool DTrackCandidate_factory::MatchMethod11(double q,DVector3 &mypos,
                   DVector3 &mymom,
                   DHelicalFit &fit2,
                   const DFDCSegment *segment1,
                   const DFDCSegment *segment2
                   ){
  if (DEBUG_LEVEL>0){
    _DBG_ << "Attempting Match method #11..." << endl;
  }

  // Package numbers
  int pack1=segment1->package;
  int pack2=segment2->package;

  // Set direction of propagation for traversing from segment1 to segment2
  if ((pack2<pack1 && mymom.z()>0.) || (pack1>pack2 && mymom.z()<0.)) mymom=(-1.)*mymom;

  // Find the magnetic field at the first hit of the second segment
  double x=segment2->hits[0]->xy.X();
  double y=segment2->hits[0]->xy.Y();
  double z=segment2->hits[0]->wire->origin.z();

  double Bz=fabs(bfield->GetBz(x,y,z));

  // Make local copies of input momentum and position
  DVector3 mymom1(mymom);
  DVector3 mypos1(mypos);
   
  // Get position and momentum for the second segment
  DVector3 mypos2,mymom2;
  GetPositionAndMomentum(z,fit2,Bz,mypos2,mymom2);

  if (pack2>pack1) mymom2=(-1.)*mymom2;

  // Swim from the first segment to the second segment using the stepper
  const DFDCPseudo *secondhit=segment2->hits[0];
  DVector3 norm(0.,0.,1.);
  stepper->SetCharge(q);
  stepper->SwimToPlane(mypos1,mymom1,secondhit->wire->origin,norm,NULL);

  //  _DBG_<< endl;
  double variance=1.;
  if (mypos1.Perp()<48.5){
    // Look for match at first hit
    double dx=mypos1.x()-secondhit->xy.X();
    double dy=mypos1.y()-secondhit->xy.Y();
    double d2=dx*dx+dy*dy;
    double prob=TMath::Prob(d2/variance,1);
    
    unsigned int num_match=(prob>0.01)?1:0;
    
    // Try to match more hits
    for (unsigned int i=1;i<segment2->hits.size();i++){
      const DFDCPseudo *hit=segment2->hits[i];
      
      if (mypos1.Perp()<48.5){
   ProjectHelixToZ(hit->wire->origin.z(),q,mymom1,mypos1);
   
   DVector2 XY=hit->xy;
   double dx=XY.X()-mypos1.x();
   double dy=XY.Y()-mypos1.y();
   double dr2=dx*dx+dy*dy;
   
   double variance=1.0;
   double prob = isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;
   if (prob>0.01) num_match++;
      }
    } 
    if (num_match>=3){ 
      if (DEBUG_LEVEL>0){
   _DBG_ << "Method 11:  found match!" << endl;
      }
      return true;
    }
  }

  // Try swimming in opposite direction     
  secondhit=segment1->hits[0];

  double q2=fit2.h*FactorForSenseOfRotation;
  stepper->SetCharge(q2);
  stepper->SwimToPlane(mypos2,mymom2,secondhit->wire->origin,norm,NULL);
  //  _DBG_ << mypos2.x() << " " << mypos2.y() << endl;
  if (mypos2.Perp()<48.5){
    // Look for match at first hit
    double dx=mypos2.x()-secondhit->xy.X();
    double dy=mypos2.y()-secondhit->xy.Y();
    double d2=dx*dx+dy*dy;
    
    double prob=TMath::Prob(d2/variance,1);
    
    unsigned int num_match=(prob>0.01)?1:0;
   
    // Try to match more hits
    for (unsigned int i=1;i<segment1->hits.size();i++){
      const DFDCPseudo *hit=segment1->hits[i];
     
      if (mypos2.Perp()<48.5){
   ProjectHelixToZ(hit->wire->origin.z(),q2,mymom2,mypos2);
   
   DVector2 XY=hit->xy;
   double dx=XY.X()-mypos2.x();
   double dy=XY.Y()-mypos2.y();
   double dr2=dx*dx+dy*dy;
      
   double variance=1.0;
   double prob = isfinite(dr2) ? TMath::Prob(dr2/variance,1):0.0;
   if (prob>0.01) num_match++;  
      }
    }
    if (num_match>=3){
      if (DEBUG_LEVEL>0){
   _DBG_ << "Method 11:  found match!" << endl;
      }
      return true;
    }
  }

  return false;
}


// Match low momentum track candidates using circle centers and radii of 
// curvature
bool DTrackCandidate_factory::MatchMethod12(DTrackCandidate *can, 
                   vector<int> &forward_matches,
                   int &num_fdc_cands_remaining){ 
  if (DEBUG_LEVEL>0){ 
    _DBG_ << "Attempting matching method #12..." <<endl; 
  }
  for (unsigned int i=0;i<fdctrackcandidates.size();i++){
    if (num_fdc_cands_remaining==0) return false;
    if (forward_matches[i]) continue;

    const DTrackCandidate *fdccan = fdctrackcandidates[i];
    if (fdccan->momentum().Mag()>0.5) continue;

    // Find how close the circle centers are to each other
    double dx=can->xc-fdccan->xc;
    double dy=can->yc-fdccan->yc;
    double dr2=dx*dx+dy*dy;
    // Require circle centers, radii and angles to agree at a reasonable level
    if (dr2<4.0 && fabs((can->rc-fdccan->rc)/can->rc)<0.5
   && fabs(can->momentum().Theta()-fdccan->momentum().Theta())<0.35 // 20 degrees
   ){  
      // Get the segments belonging to the fdc candidate
      vector<const DFDCSegment *>segments;
      fdccan->GetT(segments);
      stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing);

      // Get the hits from the input candidate
      vector<const DFDCPseudo*>myfdchits;
      can->GetT(myfdchits);
      if (myfdchits.size()>0){
   stable_sort(myfdchits.begin(),myfdchits.end(),FDCHitSortByLayerincreasing);
   unsigned int last_package
     =(myfdchits[myfdchits.size()-1]->wire->layer-1)/6;
   // look for something downstream...
   if (segments[0]->package-last_package==1){
     vector<const DCDCTrackHit *>cdchits;
     can->GetT(cdchits);
     stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing);
     // Variables for computing average Bz
     double Bz=0;
     unsigned int num_hits=0;
     
     // Instantiate the helical fitter to do the refit
     DHelicalFit fit; 
     // Add CDC axial hits
     for (unsigned int k=0;k<cdchits.size();k++){  
       if (cdchits[k]->is_stereo==false){
         double cov=0.213;  //guess
         const DVector3 origin=cdchits[k]->wire->origin;
         double x=origin.x(),y=origin.y(),z=origin.z();
         fit.AddHitXYZ(x,y,z,cov,cov,0.,true);
       }   
     }
     // Add the FDC hits and estimate Bz
     for (unsigned int k=0;k<segments.size();k++){
       for (unsigned int n=0;n<segments[k]->hits.size();n++){
          const DFDCPseudo *fdchit=segments[k]->hits[n];
          fit.AddHit(fdchit);
          //bfield->GetField(x,y,z,Bx,By,Bz);
          //Bz_avg-=Bz;
          Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),fdchit->wire->origin.z());
       }
       num_hits+=segments[k]->hits.size();
     }   
     for (unsigned int k=0;k<myfdchits.size();k++){
       const DFDCPseudo *fdchit=myfdchits[k];
       fit.AddHit(fdchit);
       //bfield->GetField(x,y,z,Bx,By,Bz);
       //Bz_avg-=Bz;
       Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),fdchit->wire->origin.z());
     }
     num_hits+=myfdchits.size();
     Bz=fabs(Bz)/double(num_hits);
     
     // Do the circle fit with all of the matched FDC hits
     if (fit.FitCircleRiemann(segments[0]->rc)==NOERROR){
       // Determine the polar angle
       double theta=fdccan->momentum().Theta();
       double theta_cdc=can->momentum().Theta();
       if (segments.size()==1 && theta_cdc<M_PI_4){
         double numcdc=double(cdchits.size());
         double numfdc=segments[0]->hits.size();
         theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
       }
       fit.tanl=tan(M_PI_2-theta);
       // Redo the line fit
       fit.FitLineRiemann();
       
       DVector3 myorigin;
       if (myfdchits.size()) myorigin.SetXYZ(myfdchits[0]->xy.X(),
                    myfdchits[0]->xy.Y(),
                    myfdchits[0]->wire->origin.z());
       if (cdchits.size()) myorigin=cdchits[0]->wire->origin;
       DVector3 pos=fdccan->position(),mom;
       UpdatePositionAndMomentum(fit,Bz,myorigin,pos,mom);

       // circle parameters
       can->rc=fit.r0;
       can->xc=fit.x0;
       can->yc=fit.y0;
       
       // Add additional fdc hits to the track as associated objects
       for (unsigned int m=0;m<segments.size();m++){
         for (unsigned int n=0;n<segments[m]->hits.size();n++){
      can->AddAssociatedObject(segments[m]->hits[n]);
         }
       }
       can->chisq=fit.chisq;
       can->Ndof=fit.ndof;
       Particle_t locPID = (FactorForSenseOfRotation*fit.h > 0.0) ? PiPlus : PiMinus;
       can->setPID(locPID);
       can->setMomentum(mom);
       can->setPosition(pos);
       
       forward_matches[i]=1;
       num_fdc_cands_remaining--;
       
       if (DEBUG_LEVEL>0){
         _DBG_ << "... Found match!" << endl;
       }
       return true;
     } // circle fit    
   } // look for packages downstream of those attached to candidate
      } // if we have fdc hits in the existing track candidate
      else{
   // Get the cdc hits belonging to the track
   vector<const DCDCTrackHit *>cdchits;
   can->GetT(cdchits);
   stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing);

   // Variables for computing average Bz
   double Bz=0;
   unsigned int num_hits=0;
   
   // Instantiate the helical fitter to do the refit
   DHelicalFit fit; 
   // Add CDC axial hits
   for (unsigned int k=0;k<cdchits.size();k++){ 
     if (cdchits[k]->is_stereo==false){
       double cov=0.213;  //guess
       const DVector3 origin=cdchits[k]->wire->origin;
       double x=origin.x(),y=origin.y(),z=origin.z();
       fit.AddHitXYZ(x,y,z,cov,cov,0.,true);
       }   
   }
   // Add the FDC hits and estimate Bz
   for (unsigned int k=0;k<segments.size();k++){
     for (unsigned int n=0;n<segments[k]->hits.size();n++){
       const DFDCPseudo *fdchit=segments[k]->hits[n];
       fit.AddHit(fdchit);
       //bfield->GetField(x,y,z,Bx,By,Bz);
          //Bz_avg-=Bz;
       Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),fdchit->wire->origin.z());
     }
     num_hits+=segments[k]->hits.size();
   }
   Bz=fabs(Bz)/double(num_hits);
   
   // Do the circle fit with all of the matched FDC hits
   if (fit.FitCircleRiemann(segments[0]->rc)==NOERROR){
     // Determine the polar angle
     double theta=fdccan->momentum().Theta();
     double theta_cdc=can->momentum().Theta();
     if (segments.size()==1 && theta_cdc<M_PI_4){
       double numcdc=double(cdchits.size());
         double numfdc=segments[0]->hits.size();
         theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
     }
     // recompute position and momentum
     fit.tanl=tan(M_PI_2-theta);
     DVector3 myorigin=cdchits[0]->wire->origin;
     DVector3 pos=fdccan->position(),mom;
     UpdatePositionAndMomentum(fit,Bz,cdchits[0]->wire->origin,pos,mom);

     // circle parameters
     can->rc=fit.r0;
     can->xc=fit.x0;
     can->yc=fit.y0;
     
     // Add additional fdc hits to the track as associated objects
     for (unsigned int m=0;m<segments.size();m++){
       for (unsigned int n=0;n<segments[m]->hits.size();n++){
         can->AddAssociatedObject(segments[m]->hits[n]);
       }
     }
      
     can->setPosition(pos);
     can->chisq=fit.chisq;
     can->Ndof=fit.ndof;
     Particle_t locPID = (FactorForSenseOfRotation*fit.h > 0.0) ? PiPlus : PiMinus;
     can->setPID(locPID);
     can->setMomentum(mom);
     
     forward_matches[i]=1;
     num_fdc_cands_remaining--;
       
     if (DEBUG_LEVEL>0){
       _DBG_ << "... Found match!" << endl;
     }
     return true;
   } // circle fit worked
      }
    } // matching criteria
  } // loop over existing track candidates

  return false;
}

// Attempts to match two single-segment track candidates by only using the 
// hit information to fit circles (i.e., not requiring the tracks to originate
// from the target region.) The idea is to recover some tracks arising from 
// secondary interactions or decays beyond the target.
bool DTrackCandidate_factory::MatchMethod13(unsigned int src_index,
                  const DTrackCandidate *srccan,
                  const DFDCSegment *segment,
                  vector<const DTrackCandidate*>&cands,
                  vector<int> &forward_matches){
  if (DEBUG_LEVEL>0){
    _DBG_ << "Attempting Match method #13..." << endl;
  }
  double q=srccan->charge();
  int pack1=segment->package;
  
  // Get hits from segment and redo fit
  DHelicalFit fit1;
  for (unsigned int n=0;n<segment->hits.size();n++){
    const DFDCPseudo *hit=segment->hits[n];
    fit1.AddHit(hit);
  }
  
  if (fit1.FitCircleRiemann(srccan->rc)!=NOERROR) return false;

  // Guess for theta and z from input candidates
  double theta=srccan->momentum().Theta();  
  fit1.tanl=tan(M_PI_2-theta);
  fit1.z_vertex=srccan->position().Z();
  
  // Redo line fit
  if (fit1.FitLineRiemann()!=NOERROR) return false;
  
  // Find the magnetic field at the first hit in the first segment
  double x=segment->hits[0]->xy.X();
  double y=segment->hits[0]->xy.Y();
  double z=segment->hits[0]->wire->origin.z();
  double Bz=fabs(bfield->GetBz(x,y,z));
  
  // Get position and momentum for the first segment
  DVector3 mypos,mymom;
  GetPositionAndMomentum(z,fit1,Bz,mypos,mymom);
  
  // Loop over the rest of the fdc track candidates, skipping those that have
  // already been used
  for (unsigned int k=src_index+1;k<forward_matches.size();k++){
    if (forward_matches[k]==0){
      const DTrackCandidate *can2=fdctrackcandidates[k];
      // Get the segment data
      vector<const DFDCSegment *>segments2;
      can2->GetT(segments2);
      
      int pack2=segments2[0]->package;
      if (abs(pack1-pack2)>0){
   // Get hits from the second segment and redo fit 
   DHelicalFit fit2;
   for (unsigned int n=0;n<segments2[0]->hits.size();n++){
     const DFDCPseudo *hit=segments2[0]->hits[n];
     fit2.AddHit(hit);
   }
   
   if (fit2.FitCircleRiemann(segments2[0]->rc)!=NOERROR) continue;

   // Guess for theta and z from input candidates
   theta=can2->momentum().Theta();  
   fit2.tanl=tan(M_PI_2-theta);
   fit2.z_vertex=srccan->position().Z();

   // Redo line fit
   if (fit2.FitLineRiemann()!=NOERROR) continue;
      
   // Try to match segments by swimming through the field
   if (MatchMethod11(q,mypos,mymom,fit2,segment,segments2[0])){
     forward_matches[k]=1;
     forward_matches[src_index]=1;

     // Create a new DTrackCandidate for output
     DTrackCandidate *can = new DTrackCandidate;
     
     // variables for finding <Bz>
     double Bz=0;
     int num_hits=0;
     
     // Add hits from first segment as associated objects
     for (unsigned int n=0;n<segment->hits.size();n++){
       const DFDCPseudo *fdchit=segment->hits[n];
       can->AddAssociatedObject(fdchit);
    
       Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(),
               fdchit->wire->origin.z());
       num_hits++;
     }
     
     // Add the hits from the second segment to the first fit object
     // and refit the circle.  Also add them as associated objects
     // to the new candidate.
     for (unsigned int n=0;n<segments2[0]->hits.size();n++){
       const DFDCPseudo *hit=segments2[0]->hits[n];
       fit1.AddHit(hit); 
       can->AddAssociatedObject(hit);
     
       Bz+=bfield->GetBz(hit->xy.X(),hit->xy.Y(),
               hit->wire->origin.z());
       num_hits++;        
     }
     Bz=fabs(Bz)/double(num_hits);
     
     // Initialize variables needed for output
     DVector3 mom=srccan->momentum();
     DVector3 pos=srccan->position();
     double q=srccan->charge(); 
     can->Ndof=srccan->Ndof;
     can->chisq=srccan->chisq;
     
     if (fit1.FitCircleRiemann(fit1.r0)==NOERROR){
       can->Ndof=fit1.ndof;
       can->chisq=fit1.chisq;
         
       // Redo line fit
       fit1.FitLineRiemann();
       
       // Guess charge from fit
       fit1.h=GetSenseOfRotation(fit1,segments2[0]->hits[0],
                  srccan->position());
       q=FactorForSenseOfRotation*fit1.h;
       
       // put z position just upstream of the first hit in z
       const DHFHit_t *myhit=fit1.GetHits()[0];
       pos.SetXYZ(myhit->x,myhit->y,myhit->z);
       GetPositionAndMomentum(myhit->z-1.,fit1,Bz,pos,mom);
     }
     // circle parameters
     can->rc=fit1.r0;
     can->xc=fit1.x0;
     can->yc=fit1.y0;

     can->setPID((q > 0.0) ? PiPlus : PiMinus);
     can->setPosition(pos);
     can->setMomentum(mom);
     
     trackcandidates.push_back(can);

     if (DEBUG_LEVEL>0){
       _DBG_ << "Match method #13 succeeded..." << endl;
     }
     return true;
   } // minimum number of matching hits
      } // different packages?
    } // already matched?
  } // loop over tracks

  return false;
}

// Routine to return momentum and position given the helical parameters and the
// z-component of the magnetic field
void
DTrackCandidate_factory::GetPositionAndMomentum(double z,const DHelicalFit &fit,
                  double Bz,DVector3 &pos,
                  DVector3 &mom) const{
  double xc=fit.x0;
  double yc=fit.y0;
  double rc=fit.r0;
  // Position
  double phi1=atan2(pos.y()-yc,pos.x()-xc);
  double q_over_rc_tanl=FactorForSenseOfRotation*fit.h/(rc*fit.tanl);
  double dphi_s=(pos.z()-z)*q_over_rc_tanl;
  double dphi1=phi1-dphi_s;// was -
  double x=xc+rc*cos(dphi1);
  double y=yc+rc*sin(dphi1);
  pos.SetXYZ(x,y,z);

  dphi1*=-1.;
  if (fit.h<0) dphi1+=M_PI;

  // Momentum 
  double pt=0.003*fabs(Bz)*rc; 
  double px=pt*sin(dphi1);
  double py=pt*cos(dphi1);
  double pz=pt*fit.tanl;
  mom.SetXYZ(px,py,pz);
}

// Use information from the start counter to try to correct the momentum and
// position for tracks seem to be pointing backwards but are probably pointing 
// forwards.  In these cases the circle projection intersects with a start 
// counter paddle but the z-position at the radius r just outside the start 
// counter barrel region is downstream of the nose, which does not make
// sense for a particle that is actually heading upstream...
bool DTrackCandidate_factory::TryToFlipDirection(vector<const DSCHit *>&schits,
                   DVector3 &mom,DVector3 &pos) const{
  if (schits.size()==0) return false;
  if (DEBUG_LEVEL>0){
    _DBG_ << "Attempting to flip direction of track..." << endl;
  }

  double zsc=sc_pos[0][1].z();
  if (pos.z()>zsc){
    unsigned int best_sc_sector_id=0;
    double dphi_min=1000.;
    double cand_phi=pos.Phi();
    for (unsigned int k=0;k<schits.size();k++){
      unsigned int sector_id=schits[k]->sector-1;
      double dphi=cand_phi-sc_pos[sector_id][0].Phi();
      if (dphi<-M_PI) dphi+=2.*M_PI;
      if (dphi>M_PI) dphi-=2.*M_PI;
      
      if (fabs(dphi)<dphi_min){
   best_sc_sector_id=sector_id;
   dphi_min=fabs(dphi);
      }
    }
    if (dphi_min<0.105){  
      // simplest approach: just change the direction -z -> +z
      mom.SetMagThetaPhi(mom.Mag(),M_PI-mom.Theta(),mom.Phi());
      pos=sc_pos[best_sc_sector_id][1];
      //     _DBG_ << " Event " << eventnumber << endl;
      if (DEBUG_LEVEL>0){
   _DBG_ << "Changing direction of CDC candidate..." << endl;
      }
      return true;
    }
  }
  

  return false;
}

// Last ditch attempt to match stray segments in the FDC to other track stubs
// after all other matching techniques have been tried...
bool DTrackCandidate_factory::MatchStraySegments(vector<int> &forward_matches,
                   int &num_fdc_cands_remaining){
  if (DEBUG_LEVEL>0){
    _DBG_ << "Attempting to match stray FDC segments..." << endl;
  }

  bool got_new_match=false;  
  for (unsigned int j=0;j<forward_matches.size();j++){
    if (num_fdc_cands_remaining==0) break;
    if (forward_matches[j]==0){
      const DTrackCandidate *srccan=fdctrackcandidates[j];
      
      // Get the segment data
      vector<const DFDCSegment *>segments;
      srccan->GetT(segments);
      const DFDCSegment *segment=segments[0];
      
      // redo circle fits for segments, forcing the circles to go through (0,0)
      if (MatchMethod9(j,srccan,segment,fdctrackcandidates,forward_matches)){
   if (DEBUG_LEVEL>0) _DBG_ << "Matched FDC segments using method #9" << endl;
   num_fdc_cands_remaining-=2;
   got_new_match=true;
      }
      // Redo circle fit assuming track is almost straight
      else if (MatchMethod10(j,srccan,segment,fdctrackcandidates,
              forward_matches)){
   if (DEBUG_LEVEL>0)  _DBG_ << "Matched FDC segments using method #10" << endl;
   num_fdc_cands_remaining-=2;
   got_new_match=true;
      }
      // Redo circle fits without any assumption that the tracks originate
      // from the target
      else if (MatchMethod13(j,srccan,segment,fdctrackcandidates,
              forward_matches)){
   if (DEBUG_LEVEL>0)  _DBG_ << "Matched FDC segments using method #13" << endl;
   num_fdc_cands_remaining-=2;
   got_new_match=true;
      }
    } // not already matched
  }
  return got_new_match;
}

// Routine for updating the position and radius given an input position pos.
// If the circle from the track intersects the circle with a radius just outside
// the start counter, returns the position and momentum at the place where the 
// two circles intersect.  If the two circles do not intersect, returns the
// position at the doca to the beam line.  If the z position is far upstream 
// of the active volume for either case, places the position at z=0.
void DTrackCandidate_factory::UpdatePositionAndMomentum(DHelicalFit &fit,
                     double Bz,
                     const DVector3 &origin,
                     DVector3 &pos,
                     DVector3 &mom) const{
  // Get position at fixed radius with respect to the beam line
  if (GetPositionAndMomentum(fit,Bz,origin,pos,mom)!=NOERROR){      
    // Get position and momentum at doca to beam line
    GetPositionAndMomentum(fit,Bz,pos,mom);
  }
  // if the z-position is far away from the active volume of the detector,
  // place position at fixed z=0.
  if (pos.z()<0){ 
    GetPositionAndMomentum(0.,fit,Bz,pos,mom);
  }
}

// Returns false if the z position at the doca to the beam line is less than 0
bool DTrackCandidate_factory::CheckZPosition(const DTrackCandidate *fdccan) 
  const{
  double phi0=atan2(-fdccan->xc,fdccan->yc);  
  double q=fdccan->charge();
  if (FactorForSenseOfRotation*q<0) phi0+=M_PI;
  double sinphi0=sin(phi0);
  double sign=(sinphi0>0)?1.:-1.;
  if (fabs(sinphi0)<1e-8) sinphi0=sign*1e-8;    
  double D=q*fdccan->rc-fdccan->xc/sinphi0;
  double x=-D*sinphi0;
  double y=D*cos(phi0);
  DVector3 pos=fdccan->position();
  double dx=pos.x()-x;
  double dy=pos.y()-y;
  double ratio=sqrt(dx*dx+dy*dy)/(2.*fdccan->rc);
  double phi_s=(ratio<1.)?2.*asin(ratio):M_PI;
  double newz=pos.z()-phi_s*tan(M_PI_2-fdccan->momentum().Theta())*fdccan->rc;

  return (newz>0);
}