/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86

Bug Summary

File:	alld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h
Warning:	line 398, column 7 Null pointer passed to 2nd parameter expecting 'nonnull'

Annotated Source Code

Press '?' to see keyboard shortcuts

Show analyzer invocation

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -analyze -disable-free -main-file-name DTrackCandidate_factory.cc -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -fno-rounding-math -mconstructor-aliases -munwind-tables -target-cpu x86-64 -tune-cpu generic -fno-split-dwarf-inlining -debugger-tuning=gdb -resource-dir /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0 -D HAVE_CCDB -D HAVE_RCDB -D HAVE_EVIO -D HAVE_TMVA=1 -D RCDB_MYSQL=1 -D RCDB_SQLITE=1 -D SQLITE_USE_LEGACY_STRUCT=ON -I .Linux_CentOS7.7-x86_64-gcc4.8.5/libraries/TRACKING -I libraries/TRACKING -I . -I libraries -I libraries/include -I /w/halld-scifs17exp/home/sdobbs/clang/halld_recon/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I external/xstream/include -I /usr/include/tirpc -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/root/root-6.08.06/include -I /w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/ccdb/ccdb_1.06.06/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/rcdb/rcdb_0.06.00/cpp/include -I /usr/include/mysql -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlitecpp/SQLiteCpp-2.2.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlite/sqlite-3.13.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/hdds/hdds-4.9.0/Linux_CentOS7.7-x86_64-gcc4.8.5/src -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/xerces-c/xerces-c-3.1.4/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/evio/evio-4.4.6/Linux-x86_64/include -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5 -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/x86_64-redhat-linux -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/backward -internal-isystem /usr/local/include -internal-isystem /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -std=c++11 -fdeprecated-macro -fdebug-compilation-dir /home/sdobbs/work/clang/halld_recon/src -ferror-limit 19 -fgnuc-version=4.2.1 -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -o /tmp/scan-build-2021-01-21-110224-160369-1 -x c++ libraries/TRACKING/DTrackCandidate_factory.cc

libraries/TRACKING/DTrackCandidate_factory.cc

→

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

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

37#include <cmath>
38using namespace std;

40#include "DTrackCandidate_factory.h"
41#include "DTrackCandidate_factory_CDC.h"
42#include "START_COUNTER/DSCHit.h"
43#include "DANA/DApplication.h"
44#include <JANA/JCalibration.h>
45#include <TRACKING/DHoughFind.h>

47#include <TROOT.h>
48#include <TH2F.h>
49#include <TMath.h>

51#define CUT10. 10.
52#define RADIUS_CUT50.0 50.0
53#define BEAM_VAR1. 1. // cm^2
54#define EPS0.001 0.001

56//------------------
57// cdc_fdc_match
58//------------------
59inline 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;
68}

70//------------------
71// SegmentSortByLayerincreasing
72//------------------
73inline 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;
82}

84//------------------
85// CDCHitSortByLayerincreasing
86//------------------
87inline 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;
99}

101//------------------
102// FDCHitSortByLayerincreasing
103//------------------
104inline 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;
115}

117//------------------
118// init
119//------------------
120jerror_t DTrackCandidate_factory::init(void)
121{
MAX_NUM_TRACK_CANDIDATES = 20;

return NOERROR;
125}

127//------------------
128// brun
129//------------------
130jerror_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.;
1
Assuming the condition is false→
2
←
'?' condition is false→

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

dIsNoFieldFlag = (dynamic_cast<const DMagneticFieldMapNoField*>(bfield) != NULL__null);
3
←
Assuming field 'bfield' is equal to NULL→
if(dIsNoFieldFlag3.1
Field 'dIsNoFieldFlag' is false
1
Field 'dIsNoFieldFlag' is false
1
Field 'dIsNoFieldFlag' is false
)
4
←
Taking false branch→
{
  //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__null;
eventLoop->GetSingle(dParticleID);
5
←
Calling 'JEventLoop::GetSingle'→

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(ROOT::GetROOT())->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(ROOT::GetROOT())->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;
209}

211//------------------
212// erun
213//------------------
214jerror_t DTrackCandidate_factory::erun(void)
215{
if (stepper) {
  delete stepper;
  stepper = nullptr;
}

return NOERROR;
222}

224//------------------
225// erun
226//------------------
227jerror_t DTrackCandidate_factory::fini(void)
228{
if (stepper) {
  delete stepper;
  stepper = nullptr;
}

return NOERROR;
235}



239//------------------
240// evnt
241//------------------
242jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
243{
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_21.57079632679489661923){
// 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__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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<360<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<377<<" " << "... matched to FDC candidate #" << i <<endl;

got_match=true;
    }
    else{  
if (DEBUG_LEVEL>0) _DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<382<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<392<<" " << "... matched to FDC candidate #" << i <<endl;
   // mark the cdc candidate as matched
   cdc_forward_matches[j]=1;
   
   if (DEBUG_LEVEL>0)
     _DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<397<<" " << ".. 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<423<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<434<<" " << "... 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_21.57079632679489661923 && !sc_pos.empty()){
if (TryToFlipDirection(schits,mom,pos)){
 if (DEBUG_LEVEL>0){
   _DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<516<<" "<< "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__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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<529<<" " << "... 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_PI3.14159265358979323846/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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<560<<" " << ".. 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<746<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<785<<" " << "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;
808}


811// Obtain position and momentum at the exit of a given package using the 
812// helical track model.
813void 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);
847}

849// Get position and momentum at doca to beam line
850void 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_PI3.14159265358979323846;
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_PI3.14159265358979323846;
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);
874}

876// Get the position and momentum at a fixed radius from the beam line
877jerror_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)>EPS0.001){
  phi_plus*=-1.;
  y_plus=yc+rc*sin(phi_plus);
}
if (fabs(r2-r2_minus)>EPS0.001){
  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_PI3.14159265358979323846-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_21.57079632679489661923;
pos.SetZ(pos0.z()-sperp*tanl);

return NOERROR;
958}

960// Find the position along a helical path at the z-position z
961void 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_21.57079632679489661923-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);
983}


986// Routine to do a crude match between cdc wires and a helical approximation to
987// the trajectory
988double 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_21.57079632679489661923-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;
1036}


1039// Find the particle charge given the helical fit result and an fdc hit 
1040// on the track
1041double 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);
1067}



1071// Redo the circle fit using all of the cdc axial wires and fdc hits associated
1072// with the track candidate.  Also compute the average Bz.
1073jerror_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;   
1117}

1119// Method to match FDC and CDC candidates that projects the CDC candidate to the
1120// cdc endplate, where the FDC candidates are reported.
1121bool 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_PI3.14159265358979323846*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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1169<<" " << ".. 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;
1180}

1182// Create new candidate after using Match Method 1 to match cdc and fdc 
1183// candidates
1184bool 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1204<<" " << "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_40.78539816339744830962){
    double numcdc=double(cdchits.size());
    double numfdc=segments[0]->hits.size();
    theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
  }
  fit.tanl=tan(M_PI_21.57079632679489661923-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;
1260}

1262// Method to match CDC and FDC candidates that projects an FDC candidate into 
1263// the CDC region using a helical approximation to the trajectory.  The code
1264// computes a doca to each wire in the cdc candidate and counts matches based
1265// 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1308<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1332<<" " << "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_40.78539816339744830962){
     double numcdc=double(cdchits.size());
     double numfdc=segments[0]->hits.size();
     theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
   }
   fit.tanl=tan(M_PI_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1446<<" " << "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_40.78539816339744830962){
 double numcdc=double(cdchits.size());
 double numfdc=segments[0]->hits.size();
 theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
      }
      fit.tanl=tan(M_PI_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1555<<" " << ".. 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;
}

1567// This method tries to match unmatched fdc candidates 
1568bool 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1578<<" " << "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__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_VAR1.,BEAM_VAR1.,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_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1724<<" " << "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;
1733}

1735// Matching method to try to deal with CDC candidates that do not point 
1736// directly at the the cdc endplate but could possibly be matched to FDC 
1737// candidates.
1738bool 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1762<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1786<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1793<<" " << "... 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_21.57079632679489661923-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;
1838}

1840// Method to match FDC candidates with stray CDC hits unassociated with track
1841// candidates
1842void 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_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1947<<" " << "... matched stray CDC hits ..." << endl;
} // match at least one cdc hit
1949}

1951// This method tries to match unmatched fdc candidates to candidates that
1952// already have fdc/cdc matches 
1953bool DTrackCandidate_factory::MatchMethod7(DTrackCandidate *srccan, 
			   vector<int> &forward_matches,
			   int &num_fdc_cands_remaining){ 
if (DEBUG_LEVEL>0)  _DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<1956<<" " << "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__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_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2104<<" " << "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;
2112}


2115// This method tries to improve the quality of the circle fit of a cdc
2116// candidate by assuming that if the outermost hit is in a stereo straw, 
2117// if this track is going to match to an fdc candidate, the location of the 
2118// hit along the straw is likely to be near the end of the straw, thus 
2119// providing another useful point on the circle.  After refitting the circle,
2120// the code adjusts the dip angle and tries to match to the remaining fdc 
2121// candidates.
2122bool 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_21.57079632679489661923);
  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__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_VAR1.,BEAM_VAR1.,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_40.78539816339744830962){
     double numcdc=double(cdchits.size());
     double numfdc=segments[0]->hits.size();
     theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
   }
   fit.tanl=tan(M_PI_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2309<<" " << "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;
2319}


2322// Method to try to match unmatched FDC segments by forcing the circle fits 
2323// to go through the origin.  
2324bool 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2330<<" " << "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_21.57079632679489661923-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_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2477<<" " << "Match method #9 succeeded..." << endl;
 }

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

return false;
2487}

2489// Method to try to match unmatched FDC segments by assuming that the tracks
2490// are sufficiently stiff we are better served by something closer to a
2491// "straight-line" approximation
2492bool 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2498<<" " << "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_21.57079632679489661923-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_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2633<<" " << "Match method #10 succeeded..." << endl;
 }
 return true;
} // minimum number of matching hits
    } // different packages?
  } // already matched?
} // loop over tracks

return false;
2642}

2644// Swims from one FDC segment to another segment looking for a match.  If this
2645// fails, swims from the second second to the first in the opposite direction. 
2646bool DTrackCandidate_factory::MatchMethod11(double q,DVector3 &mypos,
			    DVector3 &mymom,
			    DHelicalFit &fit2,
			    const DFDCSegment *segment1,
			    const DFDCSegment *segment2
			    ){
if (DEBUG_LEVEL>0){
  _DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2653<<" " << "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__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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2716<<" " << "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__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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2758<<" " << "Method 11:  found match!" << endl;
    }
    return true;
  }
}

return false;
2765}


2768// Match low momentum track candidates using circle centers and radii of 
2769// curvature
2770bool DTrackCandidate_factory::MatchMethod12(DTrackCandidate *can, 
			    vector<int> &forward_matches,
			    int &num_fdc_cands_remaining){ 
if (DEBUG_LEVEL>0){ 
  _DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2774<<" " << "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_40.78539816339744830962){
     double numcdc=double(cdchits.size());
     double numfdc=segments[0]->hits.size();
     theta=(theta*numfdc+theta_cdc*numcdc)/(numfdc+numcdc);
   }
   fit.tanl=tan(M_PI_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2888<<" " << "... 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_40.78539816339744830962){
   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_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2967<<" " << "... Found match!" << endl;
 }
 return true;
} // circle fit worked
    }
  } // matching criteria
} // loop over existing track candidates

return false;
2976}

2978// Attempts to match two single-segment track candidates by only using the 
2979// hit information to fit circles (i.e., not requiring the tracks to originate
2980// from the target region.) The idea is to recover some tracks arising from 
2981// secondary interactions or decays beyond the target.
2982bool 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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<2988<<" " << "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_21.57079632679489661923-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_21.57079632679489661923-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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<3120<<" " << "Match method #13 succeeded..." << endl;
 }
 return true;
} // minimum number of matching hits
    } // different packages?
  } // already matched?
} // loop over tracks

return false;
3129}

3131// Routine to return momentum and position given the helical parameters and the
3132// z-component of the magnetic field
3133void
3134DTrackCandidate_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_PI3.14159265358979323846;

// 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);
3158}

3160// Use information from the start counter to try to correct the momentum and
3161// position for tracks seem to be pointing backwards but are probably pointing 
3162// forwards.  In these cases the circle projection intersects with a start 
3163// counter paddle but the z-position at the radius r just outside the start 
3164// counter barrel region is downstream of the nose, which does not make
3165// sense for a particle that is actually heading upstream...
3166bool DTrackCandidate_factory::TryToFlipDirection(vector<const DSCHit *>&schits,
				 DVector3 &mom,DVector3 &pos) const{
if (schits.size()==0) return false;
if (DEBUG_LEVEL>0){
  _DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<3170<<" " << "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_PI3.14159265358979323846) dphi+=2.*M_PI3.14159265358979323846;
    if (dphi>M_PI3.14159265358979323846) dphi-=2.*M_PI3.14159265358979323846;
    
    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_PI3.14159265358979323846-mom.Theta(),mom.Phi());
    pos=sc_pos[best_sc_sector_id][1];
    //	    _DBG_ << " Event " << eventnumber << endl;
    if (DEBUG_LEVEL>0){
_DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<3195<<" " << "Changing direction of CDC candidate..." << endl;
    }
    return true;
  }
}


return false;
3203}

3205// Last ditch attempt to match stray segments in the FDC to other track stubs
3206// after all other matching techniques have been tried...
3207bool DTrackCandidate_factory::MatchStraySegments(vector<int> &forward_matches,
				 int &num_fdc_cands_remaining){
if (DEBUG_LEVEL>0){
  _DBG_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<3210<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<3226<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<3233<<" " << "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_std::cerr<<"libraries/TRACKING/DTrackCandidate_factory.cc"
<<":"<<3241<<" " << "Matched FDC segments using method #13" << endl;
num_fdc_cands_remaining-=2;
got_new_match=true;
    }
  } // not already matched
}
return got_new_match;
3248}

3250// Routine for updating the position and radius given an input position pos.
3251// If the circle from the track intersects the circle with a radius just outside
3252// the start counter, returns the position and momentum at the place where the 
3253// two circles intersect.  If the two circles do not intersect, returns the
3254// position at the doca to the beam line.  If the z position is far upstream 
3255// of the active volume for either case, places the position at z=0.
3256void 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);
}
3271}

3273// Returns false if the z position at the doca to the beam line is less than 0
3274bool DTrackCandidate_factory::CheckZPosition(const DTrackCandidate *fdccan) 
const{
double phi0=atan2(-fdccan->xc,fdccan->yc);  
double q=fdccan->charge();
if (FactorForSenseOfRotation*q<0) phi0+=M_PI3.14159265358979323846;
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_PI3.14159265358979323846;
double newz=pos.z()-phi_s*tan(M_PI_21.57079632679489661923-fdccan->momentum().Theta())*fdccan->rc;

return (newz>0);
3293}



←

/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h

→

1// $Id: JEventLoop.h 1763 2006-05-10 14:29:25Z davidl $
2//
3//    File: JEventLoop.h
4// Created: Wed Jun  8 12:30:51 EDT 2005
5// Creator: davidl (on Darwin wire129.jlab.org 7.8.0 powerpc)
6//
7 
8#ifndef _JEventLoop_
9#define _JEventLoop_
10 
11#include <sys/time.h>
12 
13#include <vector>
14#include <list>
15#include <string>
16#include <utility>
17#include <typeinfo>
18#include <string.h>
19#include <map>
20#include <utility>
21using std::vector;
22using std::list;
23using std::string;
24using std::type_info;
25 
26#include <JANA/jerror.h>
27#include <JANA/JObject.h>
28#include <JANA/JException.h>
29#include <JANA/JEvent.h>
30#include <JANA/JThread.h>
31#include <JANA/JFactory_base.h>
32#include <JANA/JCalibration.h>
33#include <JANA/JGeometry.h>
34#include <JANA/JResourceManager.h>
35#include <JANA/JStreamLog.h>
36 
37// The following is here just so we can use ROOT's THtml class to generate documentation.
38#include "cint.h"
39 
40 
41// Place everything in JANA namespace
42namespace jana{
43 
44 
45template<class T> class JFactory;
46class JApplication;
47class JEventProcessor;
48 
49 
50class JEventLoop{
51	public:
52	
53		friend class JApplication;
54	
55		enum data_source_t{
56			DATA_NOT_AVAILABLE = 1,
57			DATA_FROM_CACHE,
58			DATA_FROM_SOURCE,
59			DATA_FROM_FACTORY
60		};
61		
62		typedef struct{
63			string caller_name;
64			string caller_tag;
65			string callee_name;
66			string callee_tag;
67			double start_time;
68			double end_time;
69			data_source_t data_source;
70		}call_stack_t;
71		
72		typedef struct{
73			const char* factory_name;
74			string tag;
75			const char* filename;
76			int line;
77		}error_call_stack_t;
78 
79	                                   JEventLoop(JApplication *app); ///< Constructor
80							   virtual ~JEventLoop(); ////< Destructor
81				   virtual const char* className(void){return static_className();}
82					static const char* static_className(void){return "JEventLoop";}
83		
84						 JApplication* GetJApplication(void) const {return app;} ///< Get pointer to the JApplication object
85                                  void RefreshProcessorListFromJApplication(void); ///< Re-copy the list of JEventProcessors from JApplication
86                      virtual jerror_t AddFactory(JFactory_base* factory); ///< Add a factory
87							  jerror_t RemoveFactory(JFactory_base* factory); ///< Remove a factory
88                        JFactory_base* GetFactory(const string data_name, const char *tag="", bool allow_deftag=true); ///< Get a specific factory pointer
89                vector<JFactory_base*> GetFactories(void) const {return factories;} ///< Get all factory pointers
90                                  void GetFactoryNames(vector<string> &factorynames); ///< Get names of all factories in name:tag format
91                                  void GetFactoryNames(map<string,string> &factorynames); ///< Get names of all factories in map with key=name, value=tag
92                    map<string,string> GetDefaultTags(void) const {return default_tags;}
93                              jerror_t ClearFactories(void); ///< Reset all factories in preparation for next event.
94							  jerror_t PrintFactories(int sparsify=0); ///< Print a list of all factories.
95                              jerror_t Print(const string data_name, const char *tag=""); ///< Print the data of the given type
96 
97						 JCalibration* GetJCalibration();
98                template<class T> bool GetCalib(string namepath, map<string,T> &vals);
99                template<class T> bool GetCalib(string namepath, vector<T> &vals);
100                template<class T> bool GetCalib(string namepath, T &val);
101 
102                            JGeometry* GetJGeometry();
103                template<class T> bool GetGeom(string namepath, map<string,T> &vals);
104			    template<class T> bool GetGeom(string namepath, T &val);
105 
106                     JResourceManager* GetJResourceManager(void);
107                                string GetResource(string namepath);
108                template<class T> bool GetResource(string namepath, T vals, int event_number=0);
109 
110                                  void Initialize(void); ///< Do initializations just before event processing starts
111                              jerror_t Loop(void); ///< Loop over events
112                              jerror_t OneEvent(uint64_t event_number); ///< Process a specific single event (if source supports it)
113                              jerror_t OneEvent(void); ///< Process a single event
114                           inline void Pause(void){pause = 1;} ///< Pause event processing
115                           inline void Resume(void){pause = 0;} ///< Resume event processing
116                           inline void Quit(void){quit = 1;} ///< Clean up and exit the event loop
117                           inline bool GetQuit(void) const {return quit;}
118                                  void QuitProgram(void);
119 
120		                               // Support for random access of events
121		                          bool HasRandomAccess(void);
122		                          void AddToEventQueue(uint64_t event_number){ next_events_to_process.push_back(event_number); }
123		                          void AddToEventQueue(list<uint64_t> &event_numbers) { next_events_to_process.insert(next_events_to_process.end(), event_numbers.begin(), event_numbers.end()); }
124		                list<uint64_t> GetEventQueue(void){ return next_events_to_process; }
125		                          void ClearEventQueue(void){ next_events_to_process.clear(); }
126 
127        template<class T> JFactory<T>* GetSingle(const T* &t, const char *tag="", bool exception_if_not_one=true); ///< Get pointer to first data object from (source or factory).
128        template<class T> JFactory<T>* Get(vector<const T*> &t, const char *tag="", bool allow_deftag=true); ///< Get data object pointers from (source or factory)
129        template<class T> JFactory<T>* GetFromFactory(vector<const T*> &t, const char *tag="", data_source_t &data_source=null_data_source, bool allow_deftag=true); ///< Get data object pointers from factory
130            template<class T> jerror_t GetFromSource(vector<const T*> &t, JFactory_base *factory=NULL__null); ///< Get data object pointers from source.
131                        inline JEvent& GetJEvent(void){return event;} ///< Get pointer to the current JEvent object.
132                           inline void SetJEvent(JEvent *event){this->event = *event;} ///< Set the JEvent pointer.
133                           inline void SetAutoFree(int auto_free){this->auto_free = auto_free;} ///< Set the Auto-Free flag on/off
134                      inline pthread_t GetPThreadID(void) const {return pthread_id;} ///< Get the pthread of the thread to which this JEventLoop belongs
135                                double GetInstantaneousRate(void) const {return rate_instantaneous;} ///< Get the current event processing rate
136                                double GetIntegratedRate(void) const {return rate_integrated;} ///< Get the current event processing rate
137                                double GetLastEventProcessingTime(void) const {return delta_time_single;}
138                          unsigned int GetNevents(void) const {return Nevents;}
139 
140                           inline bool CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB);
141 
142                           inline bool GetCallStackRecordingStatus(void){ return record_call_stack; }
143                           inline void DisableCallStackRecording(void){ record_call_stack = false; }
144                           inline void EnableCallStackRecording(void){ record_call_stack = true; }
145		                     inline void CallStackStart(JEventLoop::call_stack_t &cs, const string &caller_name, const string &caller_tag, const string callee_name, const string callee_tag);
146		                     inline void CallStackEnd(JEventLoop::call_stack_t &cs);
147           inline vector<call_stack_t> GetCallStack(void){return call_stack;} ///< Get the current factory call stack
148                           inline void AddToCallStack(call_stack_t &cs){if(record_call_stack) call_stack.push_back(cs);} ///< Add specified item to call stack record but only if record_call_stack is true
149                           inline void AddToErrorCallStack(error_call_stack_t &cs){error_call_stack.push_back(cs);} ///< Add layer to the factory call stack
150     inline vector<error_call_stack_t> GetErrorCallStack(void){return error_call_stack;} ///< Get the current factory error call stack
151                                  void PrintErrorCallStack(void); ///< Print the current factory call stack
152 
153                        const JObject* FindByID(JObject::oid_t id); ///< Find a data object by its identifier.
154            template<class T> const T* FindByID(JObject::oid_t id); ///< Find a data object by its type and identifier
155                        JFactory_base* FindOwner(const JObject *t); ///< Find the factory that owns a data object by pointer
156                        JFactory_base* FindOwner(JObject::oid_t id); ///< Find a factory that owns a data object by identifier
157 
158                                       // User defined references
159                template<class T> void SetRef(T *t);    ///< Add a user reference to this JEventLoop (must be a pointer)
160                  template<class T> T* GetRef(void);   ///< Get a user-defined reference of a specific type
161          template<class T> vector<T*> GetRefsT(void); ///< Get all user-defined refrences of a specicif type
162     vector<pair<const char*, void*> > GetRefs(void){ return user_refs; }  ///< Get copy of full list of user-defined references
163                template<class T> void RemoveRef(T *t); ///< Remove user reference from list
164 
165									   // Convenience methods wrapping JEvent methods of same name
166		                      uint64_t GetStatus(void){return event.GetStatus();}
167		                          bool GetStatusBit(uint32_t bit){return event.GetStatusBit(bit);}
168		                          bool SetStatusBit(uint32_t bit, bool val=true){return event.SetStatusBit(bit, val);}
169		                          bool ClearStatusBit(uint32_t bit){return event.ClearStatusBit(bit);}
170		                          void ClearStatus(void){event.ClearStatus();}
171		                          void SetStatusBitDescription(uint32_t bit, string description){event.SetStatusBitDescription(bit, description);}
172		                        string GetStatusBitDescription(uint32_t bit){return event.GetStatusBitDescription(bit);}
173		                          void GetStatusBitDescriptions(map<uint32_t, string> &status_bit_descriptions){return event.GetStatusBitDescriptions(status_bit_descriptions);}
174                                string StatusWordToString(void);
175 
176	private:
177		JEvent event;
178		vector<JFactory_base*> factories;
179		vector<JEventProcessor*> processors;
180		vector<error_call_stack_t> error_call_stack;
181		vector<call_stack_t> call_stack;
182		JApplication *app;
183		JThread *jthread;
184		bool initialized;
185		bool print_parameters_called;
186		int pause;
187		int quit;
188		int auto_free;
189		pthread_t pthread_id;
190		map<string, string> default_tags;
191		vector<pair<string,string> > auto_activated_factories;
192		bool record_call_stack;
193		string caller_name;
194		string caller_tag;
195		vector<uint64_t> event_boundaries;
196		int32_t event_boundaries_run; ///< Run number boundaries were retrieved from (possbily 0)
197		list<uint64_t> next_events_to_process;
198		
199		uint64_t Nevents;			      ///< Total events processed (this thread)
200		uint64_t Nevents_rate;		   ///< Num. events accumulated for "instantaneous" rate
201		double delta_time_single;		///< Time spent processing last event
202		double delta_time_rate;			///< Integrated time accumulated "instantaneous" rate (partial number of events)
203		double delta_time;				///< Total time spent processing events (this thread)
204		double rate_instantaneous;		///< Latest instantaneous rate
205		double rate_integrated;			///< Rate integrated over all events
206   
207		static data_source_t null_data_source;
208 
209		vector<pair<const char*, void*> > user_refs;
210};
211 
212 
213// The following is here just so we can use ROOT's THtml class to generate documentation.
214#ifdef G__DICTIONARY
215typedef JEventLoop::call_stack_t call_stack_t;
216typedef JEventLoop::error_call_stack_t error_call_stack_t;
217#endif
218 
219#if !defined(__CINT__) && !defined(__CLING__)
220 
221//-------------
222// GetSingle
223//-------------
224template<class T>
225JFactory<T>* JEventLoop::GetSingle(const T* &t, const char *tag, bool exception_if_not_one)
226{
227	/// This is a convenience method that can be used to get a pointer to the single
228	/// object of type T from the specified factory. It simply calls the Get(vector<...>) method
229	/// and copies the first pointer into "t" (or NULL if something other than 1 object is returned).
230	/// 
231	/// This is intended to address the common situation in which there is an interest
232	/// in the event if and only if there is exactly 1 object of type T. If the event
233	/// has no objects of that type or more than 1 object of that type (for the specified
234	/// factory) then an exception of type "unsigned long" is thrown with the value
235	/// being the number of objects of type T. You can supress the exception by setting
236	/// exception_if_not_one to false. In that case, you will have to check if t==NULL to
237	/// know if the call succeeded.
238	vector<const T*> v;
239	JFactory<T> *fac = Get(v, tag);
6
←
Calling 'JEventLoop::Get'→
240 
241	if(v.size()!=1){
242		t = NULL__null;
243		if(exception_if_not_one) throw v.size();
244	}
245	
246	t = v[0];
247	
248	return fac;
249}
250 
251//-------------
252// Get
253//-------------
254template<class T> 
255JFactory<T>* JEventLoop::Get(vector<const T*> &t, const char *tag, bool allow_deftag)
256{
257	/// Retrieve or generate the array of objects of
258	/// type T for the curent event being processed
259	/// by this thread.
260	///
261	/// By default, preference is given to reading the
262	/// objects from the data source(e.g. file) before generating
263	/// them in the factory. A flag exists in the factory
264	/// however to change this so that the factory is
265	/// given preference.
266	///
267	/// Note that regardless of the setting of this flag,
268	/// the data are only either read in or generated once.
269	/// Ownership of the objects will always be with the
270	/// factory so subsequent calls will always return pointers to
271	/// the same data.
272	///
273	/// If the factory is called on to generate the data,
274	/// it is done by calling the factory's Get() method
275	/// which, in turn, calls the evnt() method. 
276	/// 
277	/// First, we just call the GetFromFactory() method.
278	/// It will make the initial decision as to whether
279	/// it should look in the source first or not. If
280	/// it returns NULL, then the factory couldn't be
281	/// found so we automatically try the file.
282	///
283	/// Note that if no factory exists to hold the objects
284	/// from the file, one can be created automatically
285	/// providing the <i>JANA:AUTOFACTORYCREATE</i>
286	/// configuration parameter is set.
287 
288	// Check if a tag was specified for this data type to use for the
289	// default.
290	const char *mytag = tag6.1
'tag' is not equal to NULL
6.1
'tag' is not equal to NULL
6.1
'tag' is not equal to NULL
==NULL__null ? "":tag; // protection against NULL tags
7
←
'?' condition is false→
291	if(strlen(mytag)==0 && allow_deftag7.1
'allow_deftag' is true
7.1
'allow_deftag' is true
7.1
'allow_deftag' is true
){
8
←
Taking true branch→
292		map<string, string>::const_iterator iter = default_tags.find(T::static_className());
293		if(iter!=default_tags.end())tag = iter->second.c_str();
9
←
Assuming the condition is true→
10
←
Taking true branch→
11
←
Value assigned to 'tag'→
294	}
295	
296	
297	// If we are trying to keep track of the call stack then we
298	// need to add a new call_stack_t object to the the list
299	// and initialize it with the start time and caller/callee
300	// info.
301	call_stack_t cs;
302 
303	// Optionally record starting info of call stack entry
304	if(record_call_stack) CallStackStart(cs, caller_name, caller_tag, T::static_className(), tag);
12
←
Assuming field 'record_call_stack' is false→
13
←
Taking false branch→
305 
306	// Get the data (or at least try to)
307	JFactory<T>* factory=NULL__null;
308	try{
309		factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
14
←
Passing value via 2nd parameter 'tag'→
15
←
Calling 'JEventLoop::GetFromFactory'→
310		if(!factory){
311			// No factory exists for this type and tag. It's possible
312			// that the source may be able to provide the objects
313			// but it will need a place to put them. We can create a
314			// dumb JFactory just to hold the data in case the source
315			// can provide the objects. Before we do though, make sure
316			// the user condones this via the presence of the
317			// "JANA:AUTOFACTORYCREATE" config parameter.
318			string p;
319			try{
320				gPARMS->GetParameter("JANA:AUTOFACTORYCREATE", p);
321			}catch(...){}
322			if(p.size()==0){
323				jout<<std::endl;
324				_DBG__std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<324<<std::endl;
325				jout<<"No factory of type \""<<T::static_className()<<"\" with tag \""<<tag<<"\" exists."<<std::endl;
326				jout<<"If you are reading objects from a file, I can auto-create a factory"<<std::endl;
327				jout<<"of the appropriate type to hold the objects, but this feature is turned"<<std::endl;
328				jout<<"off by default. To turn it on, set the \"JANA:AUTOFACTORYCREATE\""<<std::endl;
329				jout<<"configuration parameter. This can usually be done by passing the"<<std::endl;
330				jout<<"following argument to the program from the command line:"<<std::endl;
331				jout<<std::endl;
332				jout<<"   -PJANA:AUTOFACTORYCREATE=1"<<std::endl;
333				jout<<std::endl;
334				jout<<"Note that since the most commonly expected occurance of this situation."<<std::endl;
335				jout<<"is an error, the program will now throw an exception so that the factory."<<std::endl;
336				jout<<"call stack can be printed."<<std::endl;
337				jout<<std::endl;
338				throw exception();
339			}else{
340				AddFactory(new JFactory<T>(tag));
341				jout<<__FILE__"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"<<":"<<__LINE__341<<" Auto-created "<<T::static_className()<<":"<<tag<<" factory"<<std::endl;
342			
343				// Now try once more. The GetFromFactory method will call
344				// GetFromSource since it's empty.
345				factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
346			}
347		}
348	}catch(exception &e){
349		// Uh-oh, an exception was thrown. Add us to the call stack
350		// and re-throw the exception
351		error_call_stack_t ecs;
352		ecs.factory_name = T::static_className();
353		ecs.tag = tag;
354		ecs.filename = NULL__null;
355		error_call_stack.push_back(ecs);
356		throw e;
357	}
358	
359	// If recording the call stack, update the end_time field and add to stack
360	if(record_call_stack) CallStackEnd(cs);
361	
362	return factory;
363}
364 
365//-------------
366// GetFromFactory
367//-------------
368template<class T> 
369JFactory<T>* JEventLoop::GetFromFactory(vector<const T*> &t, const char *tag, data_source_t &data_source, bool allow_deftag)
370{
371	// We need to find the factory providing data type T with
372	// tag given by "tag". 
373	vector<JFactory_base*>::iterator iter=factories.begin();
374	JFactory<T> *factory = NULL__null;
375	string className(T::static_className());
376	
377	// Check if a tag was specified for this data type to use for the
378	// default.
379	const char *mytag = tag==NULL__null ? "":tag; // protection against NULL tags
16
←
Assuming 'tag' is equal to NULL→
17
←
Assuming pointer value is null→
18
←
'?' condition is true→
380	if(strlen(mytag)==0 && allow_deftag18.1
'allow_deftag' is true
18.1
'allow_deftag' is true
18.1
'allow_deftag' is true
){
19
←
Taking true branch→
381		map<string, string>::const_iterator iter = default_tags.find(className);
382		if(iter!=default_tags.end())tag = iter->second.c_str();
20
←
Assuming the condition is false→
21
←
Taking false branch→
383	}
384	
385	for(; iter!=factories.end(); iter++){
22
←
Calling 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
25
←
Returning from 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
26
←
Loop condition is true.  Entering loop body→
386		// It turns out a long standing bug in g++ makes dynamic_cast return
387		// zero improperly when used on objects created on one side of
388		// a dynamically shared object (DSO) and the cast occurs on the 
389		// other side. I saw bug reports ranging from 2001 to 2004. I saw
390		// saw it first-hand on LinuxEL4 using g++ 3.4.5. This is too bad
391		// since it is much more elegant (and safe) to use dynamic_cast.
392		// To avoid this problem which can occur with plugins, we check
393		// the name of the data classes are the same. (sigh)
394		//factory = dynamic_cast<JFactory<T> *>(*iter);
395		if(className == (*iter)->GetDataClassName())factory = (JFactory<T>*)(*iter);
27
←
Taking true branch→
396		if(factory == NULL__null)continue;
28
←
Assuming 'factory' is not equal to NULL→
29
←
Taking false branch→
397		const char *factag = factory->Tag()==NULL__null ? "":factory->Tag();
30
←
Assuming the condition is true→
31
←
'?' condition is true→
398		if(!strcmp(factag, tag)){
32
←
Null pointer passed to 2nd parameter expecting 'nonnull'
399			break;
400		}else{
401			factory=NULL__null;
402		}
403	}
404	
405	// If factory not found, just return now
406	if(!factory){
407		data_source = DATA_NOT_AVAILABLE;
408		return NULL__null;
409	}
410	
411	// OK, we found the factory. If the evnt() routine has already
412	// been called, then just call the factory's Get() routine
413	// to return a copy of the existing data
414	if(factory->evnt_was_called()){
415		factory->CopyFrom(t);
416		data_source = DATA_FROM_CACHE;
417		return factory;
418	}
419	
420	// Next option is to get the objects from the data source
421	if(factory->GetCheckSourceFirst()){
422		// If the object type/tag is found in the source, it
423		// will return NOERROR, even if there are zero instances
424		// of it. If it is not available in the source then it
425		// will return OBJECT_NOT_AVAILABLE.
426		
427		jerror_t err = GetFromSource(t, factory);
428		if(err == NOERROR){
429			// A return value of NOERROR means the source had the objects
430			// even if there were zero of them.(If the source had no
431			// information about the objects OBJECT_NOT_AVAILABLE would 
432			// have been returned.)
433			// The GetFromSource() call will eventually lead to a call to
434			// the GetObjects() method of the concrete class derived
435			// from JEventSource. That routine should copy the object
436			// pointers into the factory using the factory's CopyTo()
437			// method which also sets the evnt_called flag for the factory.
438			// Note also that the "t" vector is then filled with a call
439			// to the factory's CopyFrom() method in JEvent::GetObjects().
440			// All we need to do now is just set the factory pointers in
441			// the newly generated JObjects and return the factory pointer.
442 
443			factory->SetFactoryPointers();
444			data_source = DATA_FROM_SOURCE;
445 
446			return factory;
447		}
448	}
449		
450	// OK. It looks like we have to have the factory make this.
451	// Get pointers to data from the factory.
452	factory->Get(t);
453	factory->SetFactoryPointers();
454	data_source = DATA_FROM_FACTORY;
455	
456	return factory;
457}
458 
459//-------------
460// GetFromSource
461//-------------
462template<class T> 
463jerror_t JEventLoop::GetFromSource(vector<const T*> &t, JFactory_base *factory)
464{
465	/// This tries to get objects from the event source.
466	/// "factory" must be a valid pointer to a JFactory
467	/// object since that will take ownership of the objects
468	/// created by the source.
469	/// This should usually be called from JEventLoop::GetFromFactory
470	/// which is called from JEventLoop::Get. The latter will
471	/// create a dummy JFactory of the proper flavor and tag if
472	/// one does not already exist so if objects exist in the
473	/// file without a corresponding factory to create them, they
474	/// can still be used.
475	if(!factory)throw OBJECT_NOT_AVAILABLE;
476	
477	return event.GetObjects(t, factory);
478}
479 
480//-------------
481// CallStackStart
482//-------------
483inline void JEventLoop::CallStackStart(JEventLoop::call_stack_t &cs, const string &caller_name, const string &caller_tag, const string callee_name, const string callee_tag)
484{
485	/// This is used to fill initial info into a call_stack_t stucture
486	/// for recording the call stack. It should be matched with a call
487	/// to CallStackEnd. It is normally called from the Get() method
488	/// above, but may also be used by external actors to manipulate
489	/// the call stack (presumably for good and not evil).
490 
491	struct itimerval tmr;
492	getitimer(ITIMER_PROFITIMER_PROF, &tmr);
493 
494	cs.caller_name    = this->caller_name;
495	cs.caller_tag     = this->caller_tag;
496	this->caller_name = cs.callee_name = callee_name;
497	this->caller_tag  = cs.callee_tag  = callee_tag;
498	cs.start_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
499}
500 
501//-------------
502// CallStackEnd
503//-------------
504inline void JEventLoop::CallStackEnd(JEventLoop::call_stack_t &cs)
505{
506	/// Complete a call stack entry. This should be matched
507	/// with a previous call to CallStackStart which was
508	/// used to fill the cs structure.
509 
510	struct itimerval tmr;
511	getitimer(ITIMER_PROFITIMER_PROF, &tmr);
512	cs.end_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
513	caller_name = cs.caller_name;
514	caller_tag  = cs.caller_tag;
515	call_stack.push_back(cs);
516}
517 
518//-------------
519// CheckEventBoundary
520//-------------
521inline bool JEventLoop::CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB)
522{
523	/// Check whether the two event numbers span one or more boundaries
524	/// in the calibration/conditions database for the current run number.
525	/// Return true if they do and false if they don't. The first parameter
526	/// "event_numberA" is also checked if it lands on a boundary in which
527	/// case true is also returned. If event_numberB lands on a boundary,
528	/// but event_numberA does not, then false is returned.
529	///
530	/// This method is not expected to be called by a user. It is, however called,
531	/// everytime a JFactory's Get() method is called.
532 
533	// Make sure our copy of the boundaries is up to date
534	if(event.GetRunNumber()!=event_boundaries_run){
535		event_boundaries.clear(); // in case we can't get the JCalibration pointer
536		JCalibration *jcalib = GetJCalibration();
537		if(jcalib)jcalib->GetEventBoundaries(event_boundaries);
538		event_boundaries_run = event.GetRunNumber();
539	}
540	
541	// Loop over boundaries
542	for(unsigned int i=0; i<event_boundaries.size(); i++){
543		uint64_t eb = event_boundaries[i];
544		if((eb - event_numberA)*(eb - event_numberB) < 0.0 || eb==event_numberA){ // think about it ....
545			// events span a boundary or is on a boundary. Return true
546			return true;
547		}
548	}
549 
550	return false;
551}
552 
553//-------------
554// FindByID
555//-------------
556template<class T> 
557const T* JEventLoop::FindByID(JObject::oid_t id)
558{
559	/// This is a templated method that can be used in place
560	/// of the non-templated FindByID(oid_t) method if one knows
561	/// the class of the object with the specified id.
562	/// This method is faster than calling the non-templated
563	/// FindByID and dynamic_cast-ing the JObject since
564	/// this will only search the objects of factories that
565	/// produce the desired data type.
566	/// This method will cast the JObject pointer to one
567	/// of the specified type. To use this method,
568	/// a type is specified in the call as follows:
569	///
570	/// const DMyType *t = loop->FindByID<DMyType>(id);
571	
572	// Loop over factories looking for ones that provide
573	// specified data type.
574	for(unsigned int i=0; i<factories.size(); i++){
575		if(factories[i]->GetDataClassName() != T::static_className())continue;
576 
577		// This factory provides data of type T. Search it for
578		// the object with the specified id.
579		const JObject *my_obj = factories[i]->GetByID(id);
580		if(my_obj)return dynamic_cast<const T*>(my_obj);
581	}
582 
583	return NULL__null;
584}
585 
586//-------------
587// GetCalib (map)
588//-------------
589template<class T>
590bool JEventLoop::GetCalib(string namepath, map<string,T> &vals)
591{
592	/// Get the JCalibration object from JApplication for the run number of
593	/// the current event and call its Get() method to get the constants.
594 
595	// Note that we could do this by making "vals" a generic type T thus, combining
596	// this with the vector version below. However, doing this explicitly will make
597	// it easier for the user to understand how to call us.
598 
599	vals.clear();
600 
601	JCalibration *calib = GetJCalibration();
602	if(!calib){
603		_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<603<<" "<<"Unable to get JCalibration object for run "<<event.GetRunNumber()<<std::endl;
604		return true;
605	}
606	
607	return calib->Get(namepath, vals, event.GetEventNumber());
608}
609 
610//-------------
611// GetCalib (vector)
612//-------------
613template<class T> bool JEventLoop::GetCalib(string namepath, vector<T> &vals)
614{
615	/// Get the JCalibration object from JApplication for the run number of
616	/// the current event and call its Get() method to get the constants.
617 
618	vals.clear();
619 
620	JCalibration *calib = GetJCalibration();
621	if(!calib){
622		_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<622<<" "<<"Unable to get JCalibration object for run "<<event.GetRunNumber()<<std::endl;
623		return true;
624	}
625	
626	return calib->Get(namepath, vals, event.GetEventNumber());
627}
628 
629//-------------
630// GetCalib (single)
631//-------------
632template<class T> bool JEventLoop::GetCalib(string namepath, T &val)
633{
634	/// This is a convenience method for getting a single entry. It
635	/// simply calls the vector version and returns the first entry.
636	/// It returns true if the vector version returns true AND there
637	/// is at least one entry in the vector. No check is made for there
638	/// there being more than one entry in the vector.
639 
640	vector<T> vals;
641	bool ret = GetCalib(namepath, vals);
642	if(vals.empty()) return true;
643	val = vals[0];
644	
645	return ret;
646}
647 
648//-------------
649// GetGeom (map)
650//-------------
651template<class T>
652bool JEventLoop::GetGeom(string namepath, map<string,T> &vals)
653{
654	/// Get the JGeometry object from JApplication for the run number of
655	/// the current event and call its Get() method to get the constants.
656 
657	// Note that we could do this by making "vals" a generic type T thus, combining
658	// this with the vector version below. However, doing this explicitly will make
659	// it easier for the user to understand how to call us.
660 
661	vals.clear();
662 
663	JGeometry *geom = GetJGeometry();
664	if(!geom){
665		_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<665<<" "<<"Unable to get JGeometry object for run "<<event.GetRunNumber()<<std::endl;
666		return true;
667	}
668	
669	return geom->Get(namepath, vals);
670}
671 
672//-------------
673// GetGeom (atomic)
674//-------------
675template<class T> bool JEventLoop::GetGeom(string namepath, T &val)
676{
677	/// Get the JCalibration object from JApplication for the run number of
678	/// the current event and call its Get() method to get the constants.
679 
680	JGeometry *geom = GetJGeometry();
681	if(!geom){
682		_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<682<<" "<<"Unable to get JGeometry object for run "<<event.GetRunNumber()<<std::endl;
683		return true;
684	}
685	
686	return geom->Get(namepath, val);
687}
688 
689//-------------
690// SetRef
691//-------------
692template<class T>
693void JEventLoop::SetRef(T *t)
694{
695	pair<const char*, void*> p(typeid(T).name(), (void*)t);
696	user_refs.push_back(p);
697}
698 
699//-------------
700// GetResource
701//-------------
702template<class T> bool JEventLoop::GetResource(string namepath, T vals, int event_number)
703{
704	JResourceManager *resource_manager = GetJResourceManager();
705	if(!resource_manager){
706		string mess = string("Unable to get the JResourceManager object (namepath=\"")+namepath+"\")";
707		throw JException(mess);
708	}
709 
710	return resource_manager->Get(namepath, vals, event_number);
711}
712 
713//-------------
714// GetRef
715//-------------
716template<class T>
717T* JEventLoop::GetRef(void)
718{
719	/// Get a user-defined reference (a pointer)
720	for(unsigned int i=0; i<user_refs.size(); i++){
721		if(user_refs[i].first == typeid(T).name()) return (T*)user_refs[i].second;
722	}
723 
724	return NULL__null;
725}
726 
727//-------------
728// GetRefsT
729//-------------
730template<class T>
731vector<T*> JEventLoop::GetRefsT(void)
732{
733	vector<T*> refs;
734	for(unsigned int i=0; i<user_refs.size(); i++){
735		if(user_refs[i].first == typeid(T).name()){
736			refs.push_back((T*)user_refs[i].second);
737		}
738	}
739 
740	return refs;
741}
742 
743//-------------
744// RemoveRef
745//-------------
746template<class T>
747void JEventLoop::RemoveRef(T *t)
748{
749	vector<pair<const char*, void*> >::iterator iter;
750	for(iter=user_refs.begin(); iter!= user_refs.end(); iter++){
751		if(iter->second == (void*)t){
752			user_refs.erase(iter);
753			return;
754		}
755	}
756	_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<756<<" "<<" Attempt to remove user reference not in event loop!" << std::endl;
757}
758 
759 
760#endif //__CINT__  __CLING__
761 
762} // Close JANA namespace
763 
764 
765 
766#endif // _JEventLoop_
767 

←

/usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/bits/stl_iterator.h

1// Iterators -*- C++ -*-

3// Copyright (C) 2001-2013 Free Software Foundation, Inc.
4//
5// This file is part of the GNU ISO C++ Library.  This library is free
6// software; you can redistribute it and/or modify it under the
7// terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 3, or (at your option)
9// any later version.

11// This library is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14// GNU General Public License for more details.

16// Under Section 7 of GPL version 3, you are granted additional
17// permissions described in the GCC Runtime Library Exception, version
18// 3.1, as published by the Free Software Foundation.

20// You should have received a copy of the GNU General Public License and
21// a copy of the GCC Runtime Library Exception along with this program;
22// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23// <http://www.gnu.org/licenses/>.

25/*
*
* Copyright (c) 1994
* Hewlett-Packard Company
*
* Permission to use, copy, modify, distribute and sell this software
* and its documentation for any purpose is hereby granted without fee,
* provided that the above copyright notice appear in all copies and
* that both that copyright notice and this permission notice appear
* in supporting documentation.  Hewlett-Packard Company makes no
* representations about the suitability of this software for any
* purpose.  It is provided "as is" without express or implied warranty.
*
*
* Copyright (c) 1996-1998
* Silicon Graphics Computer Systems, Inc.
*
* Permission to use, copy, modify, distribute and sell this software
* and its documentation for any purpose is hereby granted without fee,
* provided that the above copyright notice appear in all copies and
* that both that copyright notice and this permission notice appear
* in supporting documentation.  Silicon Graphics makes no
* representations about the suitability of this software for any
* purpose.  It is provided "as is" without express or implied warranty.
*/

51/** @file bits/stl_iterator.h
*  This is an internal header file, included by other library headers.
*  Do not attempt to use it directly. @headername{iterator}
*
*  This file implements reverse_iterator, back_insert_iterator,
*  front_insert_iterator, insert_iterator, __normal_iterator, and their
*  supporting functions and overloaded operators.
*/

60#ifndef _STL_ITERATOR_H1
61#define _STL_ITERATOR_H1 1

63#include <bits/cpp_type_traits.h>
64#include <ext/type_traits.h>
65#include <bits/move.h>

67namespace std _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
68{
69_GLIBCXX_BEGIN_NAMESPACE_VERSION

/**
 * @addtogroup iterators
 * @{
 */

// 24.4.1 Reverse iterators
/**
 *  Bidirectional and random access iterators have corresponding reverse
 *  %iterator adaptors that iterate through the data structure in the
 *  opposite direction.  They have the same signatures as the corresponding
 *  iterators.  The fundamental relation between a reverse %iterator and its
 *  corresponding %iterator @c i is established by the identity:
 *  @code
 *      &*(reverse_iterator(i)) == &*(i - 1)
 *  @endcode
 *
 *  <em>This mapping is dictated by the fact that while there is always a
 *  pointer past the end of an array, there might not be a valid pointer
 *  before the beginning of an array.</em> [24.4.1]/1,2
 *
 *  Reverse iterators can be tricky and surprising at first.  Their
 *  semantics make sense, however, and the trickiness is a side effect of
 *  the requirement that the iterators must be safe.
*/
template<typename _Iterator>
  class reverse_iterator
  : public iterator<typename iterator_traits<_Iterator>::iterator_category,
      typename iterator_traits<_Iterator>::value_type,
      typename iterator_traits<_Iterator>::difference_type,
      typename iterator_traits<_Iterator>::pointer,
                    typename iterator_traits<_Iterator>::reference>
  {
  protected:
    _Iterator current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::difference_type	difference_type;
    typedef typename __traits_type::pointer		pointer;
    typedef typename __traits_type::reference		reference;

    /**
     *  The default constructor value-initializes member @p current.
     *  If it is a pointer, that means it is zero-initialized.
    */
    // _GLIBCXX_RESOLVE_LIB_DEFECTS
    // 235 No specification of default ctor for reverse_iterator
    reverse_iterator() : current() { }

    /**
     *  This %iterator will move in the opposite direction that @p x does.
    */
    explicit
    reverse_iterator(iterator_type __x) : current(__x) { }

    /**
     *  The copy constructor is normal.
    */
    reverse_iterator(const reverse_iterator& __x)
    : current(__x.current) { }

    /**
     *  A %reverse_iterator across other types can be copied if the
     *  underlying %iterator can be converted to the type of @c current.
    */
    template<typename _Iter>
      reverse_iterator(const reverse_iterator<_Iter>& __x)
: current(__x.base()) { }

    /**
     *  @return  @c current, the %iterator used for underlying work.
    */
    iterator_type
    base() const
    { return current; }

    /**
     *  @return  A reference to the value at @c --current
     *
     *  This requires that @c --current is dereferenceable.
     *
     *  @warning This implementation requires that for an iterator of the
     *           underlying iterator type, @c x, a reference obtained by
     *           @c *x remains valid after @c x has been modified or
     *           destroyed. This is a bug: http://gcc.gnu.org/PR51823
    */
    reference
    operator*() const
    {
_Iterator __tmp = current;
return *--__tmp;
    }

    /**
     *  @return  A pointer to the value at @c --current
     *
     *  This requires that @c --current is dereferenceable.
    */
    pointer
    operator->() const
    { return &(operator*()); }

    /**
     *  @return  @c *this
     *
     *  Decrements the underlying iterator.
    */
    reverse_iterator&
    operator++()
    {
--current;
return *this;
    }

    /**
     *  @return  The original value of @c *this
     *
     *  Decrements the underlying iterator.
    */
    reverse_iterator
    operator++(int)
    {
reverse_iterator __tmp = *this;
--current;
return __tmp;
    }

    /**
     *  @return  @c *this
     *
     *  Increments the underlying iterator.
    */
    reverse_iterator&
    operator--()
    {
++current;
return *this;
    }

    /**
     *  @return  A reverse_iterator with the previous value of @c *this
     *
     *  Increments the underlying iterator.
    */
    reverse_iterator
    operator--(int)
    {
reverse_iterator __tmp = *this;
++current;
return __tmp;
    }

    /**
     *  @return  A reverse_iterator that refers to @c current - @a __n
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator
    operator+(difference_type __n) const
    { return reverse_iterator(current - __n); }

    /**
     *  @return  *this
     *
     *  Moves the underlying iterator backwards @a __n steps.
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator&
    operator+=(difference_type __n)
    {
current -= __n;
return *this;
    }

    /**
     *  @return  A reverse_iterator that refers to @c current - @a __n
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator
    operator-(difference_type __n) const
    { return reverse_iterator(current + __n); }

    /**
     *  @return  *this
     *
     *  Moves the underlying iterator forwards @a __n steps.
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator&
    operator-=(difference_type __n)
    {
current += __n;
return *this;
    }

    /**
     *  @return  The value at @c current - @a __n - 1
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reference
    operator[](difference_type __n) const
    { return *(*this + __n); }
  };

//@{
/**
 *  @param  __x  A %reverse_iterator.
 *  @param  __y  A %reverse_iterator.
 *  @return  A simple bool.
 *
 *  Reverse iterators forward many operations to their underlying base()
 *  iterators.  Others are implemented in terms of one another.
 *
*/
template<typename _Iterator>
  inline bool
  operator==(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return __x.base() == __y.base(); }

template<typename _Iterator>
  inline bool
  operator<(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y.base() < __x.base(); }

template<typename _Iterator>
  inline bool
  operator!=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__x == __y); }

template<typename _Iterator>
  inline bool
  operator>(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y < __x; }

template<typename _Iterator>
  inline bool
  operator<=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__y < __x); }

template<typename _Iterator>
  inline bool
  operator>=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__x < __y); }

template<typename _Iterator>
  inline typename reverse_iterator<_Iterator>::difference_type
  operator-(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y.base() - __x.base(); }

template<typename _Iterator>
  inline reverse_iterator<_Iterator>
  operator+(typename reverse_iterator<_Iterator>::difference_type __n,
     const reverse_iterator<_Iterator>& __x)
  { return reverse_iterator<_Iterator>(__x.base() - __n); }

// _GLIBCXX_RESOLVE_LIB_DEFECTS
// DR 280. Comparison of reverse_iterator to const reverse_iterator.
template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator==(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return __x.base() == __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  { return __y.base() < __x.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator!=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__x == __y); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  { return __y < __x; }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__y < __x); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__x < __y); }

template<typename _IteratorL, typename _IteratorR>
376#if __cplusplus201103L >= 201103L
  // DR 685.
  inline auto
  operator-(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  -> decltype(__y.base() - __x.base())
382#else
  inline typename reverse_iterator<_IteratorL>::difference_type
  operator-(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
386#endif
  { return __y.base() - __x.base(); }
//@}

// 24.4.2.2.1 back_insert_iterator
/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator appends it to the container using
 *  push_back.
 *
 *  Tip:  Using the back_inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class back_insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /// The only way to create this %iterator is with a container.
    explicit
    back_insert_iterator(_Container& __x) : container(&__x) { }

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator doesn't really have a @a position in the
     *  container (you can think of the position as being permanently at
     *  the end, if you like).  Assigning a value to the %iterator will
     *  always append the value to the end of the container.
    */
427#if __cplusplus201103L < 201103L
    back_insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
container->push_back(__value);
return *this;
    }
434#else
    back_insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
container->push_back(__value);
return *this;
    }

    back_insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
container->push_back(std::move(__value));
return *this;
    }
448#endif

    /// Simply returns *this.
    back_insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    back_insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    back_insert_iterator
    operator++(int)
    { return *this; }
  };

/**
 *  @param  __x  A container of arbitrary type.
 *  @return  An instance of back_insert_iterator working on @p __x.
 *
 *  This wrapper function helps in creating back_insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container>
  inline back_insert_iterator<_Container>
  back_inserter(_Container& __x)
  { return back_insert_iterator<_Container>(__x); }

/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator prepends it to the container using
 *  push_front.
 *
 *  Tip:  Using the front_inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class front_insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /// The only way to create this %iterator is with a container.
    explicit front_insert_iterator(_Container& __x) : container(&__x) { }

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator doesn't really have a @a position in the
     *  container (you can think of the position as being permanently at
     *  the front, if you like).  Assigning a value to the %iterator will
     *  always prepend the value to the front of the container.
    */
517#if __cplusplus201103L < 201103L
    front_insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
container->push_front(__value);
return *this;
    }
524#else
    front_insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
container->push_front(__value);
return *this;
    }

    front_insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
container->push_front(std::move(__value));
return *this;
    }
538#endif

    /// Simply returns *this.
    front_insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    front_insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    front_insert_iterator
    operator++(int)
    { return *this; }
  };

/**
 *  @param  __x  A container of arbitrary type.
 *  @return  An instance of front_insert_iterator working on @p x.
 *
 *  This wrapper function helps in creating front_insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container>
  inline front_insert_iterator<_Container>
  front_inserter(_Container& __x)
  { return front_insert_iterator<_Container>(__x); }

/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator inserts it in the container at the
 *  %iterator's position, rather than overwriting the value at that
 *  position.
 *
 *  (Sequences will actually insert a @e copy of the value before the
 *  %iterator's position.)
 *
 *  Tip:  Using the inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;
    typename _Container::iterator iter;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /**
     *  The only way to create this %iterator is with a container and an
     *  initial position (a normal %iterator into the container).
    */
    insert_iterator(_Container& __x, typename _Container::iterator __i)
    : container(&__x), iter(__i) {}

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator maintains its own position in the
     *  container.  Assigning a value to the %iterator will insert the
     *  value into the container at the place before the %iterator.
     *
     *  The position is maintained such that subsequent assignments will
     *  insert values immediately after one another.  For example,
     *  @code
     *     // vector v contains A and Z
     *
     *     insert_iterator i (v, ++v.begin());
     *     i = 1;
     *     i = 2;
     *     i = 3;
     *
     *     // vector v contains A, 1, 2, 3, and Z
     *  @endcode
    */
628#if __cplusplus201103L < 201103L
    insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
iter = container->insert(iter, __value);
++iter;
return *this;
    }
636#else
    insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
iter = container->insert(iter, __value);
++iter;
return *this;
    }

    insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
iter = container->insert(iter, std::move(__value));
++iter;
return *this;
    }
652#endif

    /// Simply returns *this.
    insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    insert_iterator&
    operator++(int)
    { return *this; }
  };

/**
 *  @param __x  A container of arbitrary type.
 *  @return  An instance of insert_iterator working on @p __x.
 *
 *  This wrapper function helps in creating insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container, typename _Iterator>
  inline insert_iterator<_Container>
  inserter(_Container& __x, _Iterator __i)
  {
    return insert_iterator<_Container>(__x,
			 typename _Container::iterator(__i));
  }

// @} group iterators

691_GLIBCXX_END_NAMESPACE_VERSION
692} // namespace

694namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
695{
696_GLIBCXX_BEGIN_NAMESPACE_VERSION

// This iterator adapter is @a normal in the sense that it does not
// change the semantics of any of the operators of its iterator
// parameter.  Its primary purpose is to convert an iterator that is
// not a class, e.g. a pointer, into an iterator that is a class.
// The _Container parameter exists solely so that different containers
// using this template can instantiate different types, even if the
// _Iterator parameter is the same.
using std::iterator_traits;
using std::iterator;
template<typename _Iterator, typename _Container>
  class __normal_iterator
  {
  protected:
    _Iterator _M_current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::iterator_category iterator_category;
    typedef typename __traits_type::value_type  	value_type;
    typedef typename __traits_type::difference_type 	difference_type;
    typedef typename __traits_type::reference 	reference;
    typedef typename __traits_type::pointer   	pointer;

    _GLIBCXX_CONSTEXPRconstexpr __normal_iterator() : _M_current(_Iterator()) { }

    explicit
    __normal_iterator(const _Iterator& __i) : _M_current(__i) { }

    // Allow iterator to const_iterator conversion
    template<typename _Iter>
      __normal_iterator(const __normal_iterator<_Iter,
	  typename __enable_if<
    	       (std::__are_same<_Iter, typename _Container::pointer>::__value),
      _Container>::__type>& __i)
      : _M_current(__i.base()) { }

    // Forward iterator requirements
    reference
    operator*() const
    { return *_M_current; }

    pointer
    operator->() const
    { return _M_current; }

    __normal_iterator&
    operator++()
    {
++_M_current;
return *this;
    }

    __normal_iterator
    operator++(int)
    { return __normal_iterator(_M_current++); }

    // Bidirectional iterator requirements
    __normal_iterator&
    operator--()
    {
--_M_current;
return *this;
    }

    __normal_iterator
    operator--(int)
    { return __normal_iterator(_M_current--); }

    // Random access iterator requirements
    reference
    operator[](const difference_type& __n) const
    { return _M_current[__n]; }

    __normal_iterator&
    operator+=(const difference_type& __n)
    { _M_current += __n; return *this; }

    __normal_iterator
    operator+(const difference_type& __n) const
    { return __normal_iterator(_M_current + __n); }

    __normal_iterator&
    operator-=(const difference_type& __n)
    { _M_current -= __n; return *this; }

    __normal_iterator
    operator-(const difference_type& __n) const
    { return __normal_iterator(_M_current - __n); }

    const _Iterator&
    base() const
    { return _M_current; }
  };

// Note: In what follows, the left- and right-hand-side iterators are
// allowed to vary in types (conceptually in cv-qualification) so that
// comparison between cv-qualified and non-cv-qualified iterators be
// valid.  However, the greedy and unfriendly operators in std::rel_ops
// will make overload resolution ambiguous (when in scope) if we don't
// provide overloads whose operands are of the same type.  Can someone
// remind me what generic programming is about? -- Gaby

// Forward iterator requirements
template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator==(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() == __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator==(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() == __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator!=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() != __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator!=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() != __rhs.base(); }
23
←
Assuming the condition is true→
24
←
Returning the value 1, which participates in a condition later→

// Random access iterator requirements
template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator<(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() < __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator<(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() < __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator>(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() > __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator>(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() > __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator<=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() <= __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator<=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() <= __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator>=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() >= __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator>=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() >= __rhs.base(); }

// _GLIBCXX_RESOLVE_LIB_DEFECTS
// According to the resolution of DR179 not only the various comparison
// operators but also operator- must accept mixed iterator/const_iterator
// parameters.
template<typename _IteratorL, typename _IteratorR, typename _Container>
881#if __cplusplus201103L >= 201103L
  // DR 685.
  inline auto
  operator-(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  -> decltype(__lhs.base() - __rhs.base())
887#else
  inline typename __normal_iterator<_IteratorL, _Container>::difference_type
  operator-(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
891#endif
  { return __lhs.base() - __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline typename __normal_iterator<_Iterator, _Container>::difference_type
  operator-(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() - __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline __normal_iterator<_Iterator, _Container>
  operator+(typename __normal_iterator<_Iterator, _Container>::difference_type
     __n, const __normal_iterator<_Iterator, _Container>& __i)
  { return __normal_iterator<_Iterator, _Container>(__i.base() + __n); }

906_GLIBCXX_END_NAMESPACE_VERSION
907} // namespace

909#if __cplusplus201103L >= 201103L

911namespace std _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
912{
913_GLIBCXX_BEGIN_NAMESPACE_VERSION

/**
 * @addtogroup iterators
 * @{
 */

// 24.4.3  Move iterators
/**
 *  Class template move_iterator is an iterator adapter with the same
 *  behavior as the underlying iterator except that its dereference
 *  operator implicitly converts the value returned by the underlying
 *  iterator's dereference operator to an rvalue reference.  Some
 *  generic algorithms can be called with move iterators to replace
 *  copying with moving.
 */
template<typename _Iterator>
  class move_iterator
  {
  protected:
    _Iterator _M_current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::iterator_category iterator_category;
    typedef typename __traits_type::value_type  	value_type;
    typedef typename __traits_type::difference_type	difference_type;
    // NB: DR 680.
    typedef _Iterator					pointer;
    typedef value_type&&				reference;

    move_iterator()
    : _M_current() { }

    explicit
    move_iterator(iterator_type __i)
    : _M_current(__i) { }

    template<typename _Iter>
move_iterator(const move_iterator<_Iter>& __i)
: _M_current(__i.base()) { }

    iterator_type
    base() const
    { return _M_current; }

    reference
    operator*() const
    { return std::move(*_M_current); }

    pointer
    operator->() const
    { return _M_current; }

    move_iterator&
    operator++()
    {
++_M_current;
return *this;
    }

    move_iterator
    operator++(int)
    {
move_iterator __tmp = *this;
++_M_current;
return __tmp;
    }

    move_iterator&
    operator--()
    {
--_M_current;
return *this;
    }

    move_iterator
    operator--(int)
    {
move_iterator __tmp = *this;
--_M_current;
return __tmp;
    }

    move_iterator
    operator+(difference_type __n) const
    { return move_iterator(_M_current + __n); }

    move_iterator&
    operator+=(difference_type __n)
    {
_M_current += __n;
return *this;
    }

    move_iterator
    operator-(difference_type __n) const
    { return move_iterator(_M_current - __n); }
  
    move_iterator&
    operator-=(difference_type __n)
    { 
_M_current -= __n;
return *this;
    }

    reference
    operator[](difference_type __n) const
    { return std::move(_M_current[__n]); }
  };

// Note: See __normal_iterator operators note from Gaby to understand
// why there are always 2 versions for most of the move_iterator
// operators.
template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator==(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return __x.base() == __y.base(); }

template<typename _Iterator>
  inline bool
  operator==(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return __x.base() == __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator!=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__x == __y); }

template<typename _Iterator>
  inline bool
  operator!=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__x == __y); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  { return __x.base() < __y.base(); }

template<typename _Iterator>
  inline bool
  operator<(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  { return __x.base() < __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__y < __x); }

template<typename _Iterator>
  inline bool
  operator<=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__y < __x); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  { return __y < __x; }

template<typename _Iterator>
  inline bool
  operator>(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  { return __y < __x; }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__x < __y); }

template<typename _Iterator>
  inline bool
  operator>=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__x < __y); }

// DR 685.
template<typename _IteratorL, typename _IteratorR>
  inline auto
  operator-(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  -> decltype(__x.base() - __y.base())
  { return __x.base() - __y.base(); }

template<typename _Iterator>
  inline auto
  operator-(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  -> decltype(__x.base() - __y.base())
  { return __x.base() - __y.base(); }

template<typename _Iterator>
  inline move_iterator<_Iterator>
  operator+(typename move_iterator<_Iterator>::difference_type __n,
     const move_iterator<_Iterator>& __x)
  { return __x + __n; }

template<typename _Iterator>
  inline move_iterator<_Iterator>
  make_move_iterator(_Iterator __i)
  { return move_iterator<_Iterator>(__i); }

template<typename _Iterator, typename _ReturnType
  = typename conditional<__move_if_noexcept_cond
    <typename iterator_traits<_Iterator>::value_type>::value,
              _Iterator, move_iterator<_Iterator>>::type>
  inline _ReturnType
  __make_move_if_noexcept_iterator(_Iterator __i)
  { return _ReturnType(__i); }

// @} group iterators

1137_GLIBCXX_END_NAMESPACE_VERSION
1138} // namespace

1140#define _GLIBCXX_MAKE_MOVE_ITERATOR(_Iter)std::make_move_iterator(_Iter) std::make_move_iterator(_Iter)
1141#define _GLIBCXX_MAKE_MOVE_IF_NOEXCEPT_ITERATOR(_Iter)std::__make_move_if_noexcept_iterator(_Iter) \
std::__make_move_if_noexcept_iterator(_Iter)
1143#else
1144#define _GLIBCXX_MAKE_MOVE_ITERATOR(_Iter)std::make_move_iterator(_Iter) (_Iter)
1145#define _GLIBCXX_MAKE_MOVE_IF_NOEXCEPT_ITERATOR(_Iter)std::__make_move_if_noexcept_iterator(_Iter) (_Iter)
1146#endif // C++11

1148#endif