/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 DTrackFitterStraightTrack.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/DTrackFitterStraightTrack.cc

libraries/TRACKING/DTrackFitterStraightTrack.cc

→

1// $Id$
2//
3//    File: DTrackFitterStraightTrack.cc
4// Created: Tue Mar 12 10:22:32 EDT 2019
5// Creator: staylor (on Linux ifarm1402.jlab.org 3.10.0-327.el7.x86_64 x86_64)
6//

8#include "DTrackFitterStraightTrack.h"

10//---------------------------------
11// DTrackFitterStraightTrack    (Constructor)
12//---------------------------------
13DTrackFitterStraightTrack::DTrackFitterStraightTrack(JEventLoop *loop):DTrackFitter(loop){
int runnumber=(loop->GetJEvent()).GetRunNumber();
DApplication* dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
JCalibration *jcalib = dapp->GetJCalibration(runnumber);
const DGeometry *geom = dapp->GetDGeometry(runnumber);

// Get pointer to TrackFinder object 
vector<const DTrackFinder *> finders;
loop->Get(finders);
1
Calling 'JEventLoop::Get'→
  
 // Drop the const qualifier from the DTrackFinder pointer
finder = const_cast<DTrackFinder*>(finders[0]);

// Resource pool for matrices
dResourcePool_TMatrixFSym = std::make_shared<DResourcePool<TMatrixFSym>>
  ();
dResourcePool_TMatrixFSym->Set_ControlParams(20, 20, 50);

//****************
// FDC parameters
//****************

map<string,double>drift_res_parms;
jcalib->Get("FDC/drift_resolution_parms",drift_res_parms); 
DRIFT_RES_PARMS[0]=drift_res_parms["p0"];   
DRIFT_RES_PARMS[1]=drift_res_parms["p1"];
DRIFT_RES_PARMS[2]=drift_res_parms["p2"]; 

// Time-to-distance function parameters for FDC
map<string,double>drift_func_parms;
jcalib->Get("FDC/drift_function_parms",drift_func_parms); 
DRIFT_FUNC_PARMS[0]=drift_func_parms["p0"];   
DRIFT_FUNC_PARMS[1]=drift_func_parms["p1"];
DRIFT_FUNC_PARMS[2]=drift_func_parms["p2"]; 
DRIFT_FUNC_PARMS[3]=drift_func_parms["p3"];
DRIFT_FUNC_PARMS[4]=1000.;
DRIFT_FUNC_PARMS[5]=0.;
map<string,double>drift_func_ext;
if (jcalib->Get("FDC/drift_function_ext",drift_func_ext)==false){
  DRIFT_FUNC_PARMS[4]=drift_func_ext["p4"]; 
  DRIFT_FUNC_PARMS[5]=drift_func_ext["p5"]; 
}

PLANE_TO_SKIP=0;
gPARMS->SetDefaultParameter("TRKFIT:PLANE_TO_SKIP",PLANE_TO_SKIP); 

//****************
// CDC parameters
//****************
RING_TO_SKIP=0;
gPARMS->SetDefaultParameter("TRKFIT:RING_TO_SKIP",RING_TO_SKIP);

map<string, double> cdc_res_parms;
jcalib->Get("CDC/cdc_resolution_parms_v2", cdc_res_parms);
CDC_RES_PAR1 = cdc_res_parms["res_par1"];
CDC_RES_PAR2 = cdc_res_parms["res_par2"];
CDC_RES_PAR3 = cdc_res_parms["res_par3"];

vector< map<string, double> > tvals;
cdc_drift_table.clear();  
if (jcalib->Get("CDC/cdc_drift_table", tvals)==false){    
  for(unsigned int i=0; i<tvals.size(); i++){
    map<string, double> &row = tvals[i];
    cdc_drift_table.push_back(1000.*row["t"]);
  }
}
else{
  jerr << " CDC time-to-distance table not available... bailing..." << endl;
  exit(0);
}

if (jcalib->Get("CDC/drift_parameters", tvals)==false){
  map<string, double> &row = tvals[0]; //long drift side
  long_drift_func[0][0]=row["a1"]; 
  long_drift_func[0][1]=row["a2"];
  long_drift_func[0][2]=row["a3"];  
  long_drift_func[1][0]=row["b1"];
  long_drift_func[1][1]=row["b2"];
  long_drift_func[1][2]=row["b3"];
  long_drift_func[2][0]=row["c1"];
  long_drift_func[2][1]=row["c2"];
  long_drift_func[2][2]=row["c3"];
  
  row = tvals[1]; // short drift side
  short_drift_func[0][0]=row["a1"];
  short_drift_func[0][1]=row["a2"];
  short_drift_func[0][2]=row["a3"];  
  short_drift_func[1][0]=row["b1"];
  short_drift_func[1][1]=row["b2"];
  short_drift_func[1][2]=row["b3"];
  short_drift_func[2][0]=row["c1"];
  short_drift_func[2][1]=row["c2"];
  short_drift_func[2][2]=row["c3"];
}

 // Get the straw sag parameters from the database
 max_sag.clear();
 sag_phi_offset.clear();
 unsigned int numstraws[28]={42,42,54,54,66,66,80,80,93,93,106,106,123,123,
    135,135,146,146,158,158,170,170,182,182,197,197,209,209};
 unsigned int straw_count=0,ring_count=0;
 if (jcalib->Get("CDC/sag_parameters", tvals)==false){
   vector<double>temp,temp2;
   for(unsigned int i=0; i<tvals.size(); i++){
     map<string, double> &row = tvals[i];
     
     temp.push_back(row["offset"]);
     temp2.push_back(row["phi"]);
     
     straw_count++;
     if (straw_count==numstraws[ring_count]){
max_sag.push_back(temp);
sag_phi_offset.push_back(temp2);
temp.clear();
temp2.clear();
straw_count=0;
ring_count++;
     }
   }
 }

 // CDC Geometry
 vector<double>cdc_origin;
 vector<double>cdc_center;
 vector<double>cdc_upstream_endplate_pos; 
 vector<double>cdc_endplate_dim;
 geom->Get("//posXYZ[@volume='CentralDC'/@X_Y_Z",cdc_origin);
 geom->Get("//posXYZ[@volume='centralDC']/@X_Y_Z",cdc_center);
 geom->Get("//posXYZ[@volume='CDPU']/@X_Y_Z",cdc_upstream_endplate_pos);
 geom->Get("//tubs[@name='CDPU']/@Rio_Z",cdc_endplate_dim);
 cdc_origin[2]+=cdc_center[2]+cdc_upstream_endplate_pos[2]
    +0.5*cdc_endplate_dim[2];
 upstreamEndplate=cdc_origin[2];

 double endplate_dz=0.;
 geom->GetCDCEndplate(downstreamEndplate,endplate_dz,cdc_endplate_rmin,
	cdc_endplate_rmax);
 downstreamEndplate-=0.5*endplate_dz;

 // Outer detector geometry parameters
 if (geom->GetDIRCZ(dDIRCz)==false) dDIRCz=1000.;
 geom->GetFCALZ(dFCALz); 
 vector<double>tof_face;
 geom->Get("//section/composition/posXYZ[@volume='ForwardTOF']/@X_Y_Z",
    tof_face);
 vector<double>tof_plane;  
 geom->Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='0']", tof_plane);
 dTOFz=tof_face[2]+tof_plane[2]; 
 geom->Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='1']", tof_plane);
 dTOFz+=tof_face[2]+tof_plane[2];
 dTOFz*=0.5;  // mid plane between tof planes
 geom->GetTRDZ(dTRDz_vec); // TRD planes
 
 // Get start counter geometry;
 if (geom->GetStartCounterGeom(sc_pos,sc_norm)){
   // Create vector of direction vectors in scintillator planes
   for (int i=0;i<30;i++){
     vector<DVector3>temp;
     for (unsigned int j=0;j<sc_pos[i].size()-1;j++){
double dx=sc_pos[i][j+1].x()-sc_pos[i][j].x();
double dy=sc_pos[i][j+1].y()-sc_pos[i][j].y();
double dz=sc_pos[i][j+1].z()-sc_pos[i][j].z();
temp.push_back(DVector3(dx/dz,dy/dz,1.));
     }
     sc_dir.push_back(temp);
   }
   SC_END_NOSE_Z=sc_pos[0][12].z();
   SC_BARREL_R=sc_pos[0][0].Perp();
   SC_PHI_SECTOR1=sc_pos[0][0].Phi();
 }
 

 //*************************
 // Command-line parameters
 //*************************
 VERBOSE=0;
 gPARMS->SetDefaultParameter("TRKFIT:VERBOSE",VERBOSE);
 COSMICS=false;
 gPARMS->SetDefaultParameter("TRKFIND:COSMICS",COSMICS);
 CHI2CUT = 15.0; 
 gPARMS->SetDefaultParameter("TRKFIT:CHI2CUT",CHI2CUT);
 DO_PRUNING = true;
 gPARMS->SetDefaultParameter("TRKFIT:DO_PRUNING",DO_PRUNING);
196}

198//---------------------------------
199// ~DTrackFitterStraightTrack    (Destructor)
200//---------------------------------
201DTrackFitterStraightTrack::~DTrackFitterStraightTrack()
202{

204}

206//----------------------------------------------------
207// Comparison routines for sorting
208//----------------------------------------------------
209bool DTrackFitterStraightTrack_cdc_hit_cmp(const DCDCTrackHit *a,
			   const DCDCTrackHit *b){

 return(a->wire->origin.Y()>b->wire->origin.Y());
213}

215bool DTrackFitterStraightTrack_cdc_hit_reverse_cmp(const DCDCTrackHit *a,
				   const DCDCTrackHit *b){

 return(a->wire->origin.Y()<b->wire->origin.Y());
219}

221bool DTrackFitterStraightTrack_cdc_hit_radius_cmp(const DCDCTrackHit *a,
				  const DCDCTrackHit *b){
 if (a==NULL__null || b==NULL__null){
    cout << "Null pointer in CDC hit list??" << endl;
    return false;
 }
 const DCDCWire *wire_a= a->wire;
 const DCDCWire *wire_b= b->wire;
 if(wire_a->ring == wire_b->ring){
    return wire_a->straw < wire_b->straw;
 }

 return (wire_a->ring<wire_b->ring);
234}

236bool DTrackFitterStraightTrack_fdc_hit_cmp(const DFDCPseudo *a,
			   const DFDCPseudo *b){

 return(a->wire->origin.z()<b->wire->origin.z());
240}

242//----------------------------------------------------
243// CDC fitting routines
244//----------------------------------------------------

246// Smearing function derived from fitting residuals
247inline double DTrackFitterStraightTrack::CDCDriftVariance(double t) const { 
 if (t<0.) t=0.;
 double sigma=CDC_RES_PAR1/(t+1.)+CDC_RES_PAR2 + CDC_RES_PAR3*t;
 return sigma*sigma;
251}

253// Convert time to distance for the cdc.  Also returns the variance in d.
254void DTrackFitterStraightTrack::CDCDriftParameters(double dphi,double delta,
				   double t, double &d, 
				   double &V) const{

258// NSJ 26 May 2020 included long side a3, b3 and short side c1, c2, c3
259// Previously these parameters were not used (0 in ccdb) for production runs
260// except intensity scan run 72312 by accident 5 May 2020, superseded 8 May.
261// They were used in 2015 for runs 0-3050.

263d=0.; 
double dd_dt=0;
if (t>0){
  double f_0=0.;
  double f_delta=0.;

  if (delta>0){
    double a1=long_drift_func[0][0];
    double a2=long_drift_func[0][1];
    double a3=long_drift_func[0][2];
    double b1=long_drift_func[1][0];
    double b2=long_drift_func[1][1];
    double b3=long_drift_func[1][2];
    double c1=long_drift_func[2][0];
    double c2=long_drift_func[2][1];
    double c3=long_drift_func[2][2];
    
    // use "long side" functional form
    double my_t=0.001*t;
    double sqrt_t=sqrt(my_t);
    double t3=my_t*my_t*my_t;
    double delta_mag=fabs(delta);

    double delta_sq=delta*delta;
    double a=a1+a2*delta_mag+a3*delta_sq;
    double b=b1+b2*delta_mag+b3*delta_sq;
    double c=c1+c2*delta_mag+c3*delta_sq;

    f_delta=a*sqrt_t+b*my_t+c*t3;
    f_0=a1*sqrt_t+b1*my_t+c1*t3;

    dd_dt=0.001*(0.5*a/sqrt_t+b+3.*c*my_t*my_t);

    }
  else{
    double my_t=0.001*t;
    double sqrt_t=sqrt(my_t);
    double t3=my_t*my_t*my_t;
    double delta_mag=fabs(delta);
    
    // use "short side" functional form
    double a1=short_drift_func[0][0];
    double a2=short_drift_func[0][1];
    double a3=short_drift_func[0][2];
    double b1=short_drift_func[1][0];
    double b2=short_drift_func[1][1];
    double b3=short_drift_func[1][2];
    double c1=short_drift_func[2][0];
    double c2=short_drift_func[2][1];
    double c3=short_drift_func[2][2];

    double delta_sq=delta*delta;
    double a=a1+a2*delta_mag+a3*delta_sq;
    double b=b1+b2*delta_mag+b3*delta_sq;
    double c=c1+c2*delta_mag+c3*delta_sq;

    f_delta=a*sqrt_t+b*my_t+c*t3;
    f_0=a1*sqrt_t+b1*my_t+c1*t3;

    dd_dt=0.001*(0.5*a/sqrt_t+b+3.*c*my_t*my_t);
    
  }
  
  unsigned int max_index=cdc_drift_table.size()-1;
  if (t>cdc_drift_table[max_index]){
    //_DBG_ << "t: " << t <<" d " << f_delta <<endl;
    d=f_delta;
  }
  
  // Drift time is within range of table -- interpolate...
  unsigned int index=0;
  index=Locate(cdc_drift_table,t);
  double dt=cdc_drift_table[index+1]-cdc_drift_table[index];
  double frac=(t-cdc_drift_table[index])/dt;
  double d_0=0.01*(double(index)+frac); 

  double P=0.;
  double tcut=250.0; // ns
  if (t<tcut) {
    P=(tcut-t)/tcut;
  }
  d=f_delta*(d_0/f_0*P+1.-P);
}
double VarMs=0.001; // kludge for material effects
V=CDCDriftVariance(t)+mVarT0*dd_dt*dd_dt+VarMs;
348}

350// Perform the Kalman Filter for the current set of cdc hits
351jerror_t DTrackFitterStraightTrack::KalmanFilter(DMatrix4x1 &S,DMatrix4x4 &C,
				 vector<int>&used_hits,
				 vector<cdc_update_t>&updates,
				 double &chi2,
				 int &ndof,
				 unsigned int iter){
DMatrix1x4 H;  // Track projection matrix
DMatrix4x1 H_T; // Transpose of track projection matrix 
DMatrix4x1 K;  // Kalman gain matrix
DMatrix4x4 J; // Jacobian matrix
DMatrix4x1 S0; // State vector from reference trajectory
double V=1.15*(0.78*0.78/12.); // sigma=cell_size/sqrt(12.)*scale_factor

const double d_EPS=1e-8;

// Zero out the vector of used hit flags
for (unsigned int i=0;i<used_hits.size();i++) used_hits[i]=0;

//Initialize chi2 and ndof
chi2=0.;
ndof=0;

//Save the starting values for C and S
trajectory[0].Skk=S;
trajectory[0].Ckk=C;

double doca2=0.;

// CDC index and wire position variables
unsigned int cdc_index=cdchits.size()-1;
bool more_hits=true;
const DCDCWire *wire=cdchits[cdc_index]->wire;
DVector3 origin=wire->origin;
double z0=origin.z();
double vz=wire->udir.z();
if (VERBOSE) jout << " Starting in Ring " << wire->ring << endl;
DVector3 wdir=(1./vz)*wire->udir;
DVector3 wirepos=origin+(trajectory[0].z-z0)*wdir;

/// compute initial doca^2 to first wire
double dx=S(state_x)-wirepos.X();
double dy=S(state_y)-wirepos.Y();
double old_doca2=dx*dx+dy*dy;

// Loop over all steps in the trajectory
S0=trajectory[0].S;
J=trajectory[0].J;
for (unsigned int k=1;k<trajectory.size();k++){
  if (!C.IsPosDef()) return UNRECOVERABLE_ERROR;
  
  // Propagate the state and covariance matrix forward in z
  S=trajectory[k].S+J*(S-S0);
  C=J*C*J.Transpose();
  
  // Save the current state and covariance matrix in the deque
  trajectory[k].Skk=S;
  trajectory[k].Ckk=C;
  
  // Save S and J for the next step
  S0=trajectory[k].S;
  J=trajectory[k].J;
  
  // Position along wire
  wirepos=origin+(trajectory[k].z-z0)*wdir;
  
  // New doca^2
  dx=S(state_x)-wirepos.X();
  dy=S(state_y)-wirepos.Y();
  doca2=dx*dx+dy*dy;
  if (VERBOSE > 10) jout<< "At Position " << S(state_x) << " " << S(state_y) << " " << trajectory[k].z << " doca2 " << doca2 << endl;
  
  if (doca2>old_doca2 && more_hits){
     
    // zero-position and direction of line describing particle trajectory
    double tx=S(state_tx),ty=S(state_ty);
    DVector3 pos0(S(state_x),S(state_y),trajectory[k].z);
    DVector3 tdir(tx,ty,1.);
    
    // Find the true doca to the wire
    DVector3 diff=pos0-origin;
    double dx0=diff.x(),dy0=diff.y();
    double wdir_dot_diff=diff.Dot(wdir);
    double tdir_dot_diff=diff.Dot(tdir);
    double tdir_dot_wdir=tdir.Dot(wdir);
    double tdir2=tdir.Mag2();
    double wdir2=wdir.Mag2();
    double D=tdir2*wdir2-tdir_dot_wdir*tdir_dot_wdir;
    double N=tdir_dot_wdir*wdir_dot_diff-wdir2*tdir_dot_diff;
    double N1=tdir2*wdir_dot_diff-tdir_dot_wdir*tdir_dot_diff;
    double scale=1./D;
    double s=scale*N;
    double t=scale*N1;
    diff+=s*tdir-t*wdir;
    double d=diff.Mag()+d_EPS; // prevent division by zero
     
    // The next measurement
    double dmeas=0.39,delta=0.;
    double tdrift=cdchits[cdc_index]->tdrift-trajectory[k].t;
    if (fit_type==kTimeBased){
double phi_d=diff.Phi();
double dphi=phi_d-origin.Phi();  
while (dphi>M_PI3.14159265358979323846) dphi-=2*M_PI3.14159265358979323846;
while (dphi<-M_PI3.14159265358979323846) dphi+=2*M_PI3.14159265358979323846;

int ring_index=cdchits[cdc_index]->wire->ring-1;
int straw_index=cdchits[cdc_index]->wire->straw-1;
double dz=t*wdir.z();
double delta=max_sag[ring_index][straw_index]*(1.-dz*dz/5625.)
 *cos(phi_d+sag_phi_offset[ring_index][straw_index]);
CDCDriftParameters(dphi,delta,tdrift,dmeas,V);
    }

    // residual
    double res=dmeas-d;
    if (VERBOSE>5) jout << " Residual " << res << endl;
    
    // Track projection
    
    double one_over_d=1./d;
    double diffx=diff.x(),diffy=diff.y(),diffz=diff.z();
    double wx=wdir.x(),wy=wdir.y();
    
    double dN1dtx=2.*tx*wdir_dot_diff-wx*tdir_dot_diff-tdir_dot_wdir*dx0;
    double dDdtx=2.*tx*wdir2-2.*tdir_dot_wdir*wx;
    double dtdtx=scale*(dN1dtx-t*dDdtx);
    
    double dN1dty=2.*ty*wdir_dot_diff-wy*tdir_dot_diff-tdir_dot_wdir*dy0;
    double dDdty=2.*ty*wdir2-2.*tdir_dot_wdir*wy;
    double dtdty=scale*(dN1dty-t*dDdty);
    
    double dNdtx=wx*wdir_dot_diff-wdir2*dx0;
    double dsdtx=scale*(dNdtx-s*dDdtx);
    
    double dNdty=wy*wdir_dot_diff-wdir2*dy0;
    double dsdty=scale*(dNdty-s*dDdty);
    
    H(state_tx)=H_T(state_tx)
=one_over_d*(diffx*(s+tx*dsdtx-wx*dtdtx)+diffy*(ty*dsdtx-wy*dtdtx)
     +diffz*(dsdtx-dtdtx));
    H(state_ty)=H_T(state_ty)
=one_over_d*(diffx*(tx*dsdty-wx*dtdty)+diffy*(s+ty*dsdty-wy*dtdty)
     +diffz*(dsdty-dtdty));
    
    double dsdx=scale*(tdir_dot_wdir*wx-wdir2*tx);
    double dtdx=scale*(tdir2*wx-tdir_dot_wdir*tx);
    double dsdy=scale*(tdir_dot_wdir*wy-wdir2*ty);
    double dtdy=scale*(tdir2*wy-tdir_dot_wdir*ty);
    
    H(state_x)=H_T(state_x)
=one_over_d*(diffx*(1.+dsdx*tx-dtdx*wx)+diffy*(dsdx*ty-dtdx*wy)
     +diffz*(dsdx-dtdx));
    H(state_y)=H_T(state_y)
=one_over_d*(diffx*(dsdy*tx-dtdy*wx)+diffy*(1.+dsdy*ty-dtdy*wy)
     +diffz*(dsdy-dtdy));

    double InvV=1./(V+H*C*H_T);
    
    // Check how far this hit is from the projection
    double chi2check=res*res*InvV;
    if (chi2check < CHI2CUT || DO_PRUNING == 0 || (COSMICS && iter == 0)){
if (VERBOSE>5) jout << "Hit Added to track " << endl;
if (cdchits[cdc_index]->wire->ring!=RING_TO_SKIP){
 
 // Compute Kalman gain matrix
 K=InvV*(C*H_T);
 
 // Update state vector covariance matrix
 DMatrix4x4 Ctest=C-K*(H*C);
 
 //C.Print();
 //K.Print();
 //Ctest.Print();
 
 // Check that Ctest is positive definite
 if (!Ctest.IsPosDef()) return VALUE_OUT_OF_RANGE;
 C=Ctest;
 if(VERBOSE>10) C.Print();
 // Update the state vector 
 //S=S+res*K;
 S+=res*K;
 if(VERBOSE>10) S.Print();
 
 // Compute new residual 
 //d=finder->FindDoca(trajectory[k].z,S,wdir,origin);
 res=res-H*K*res;
 
 // Update chi2 
 double fit_V=V-H*C*H_T;
 chi2+=res*res/fit_V;
 ndof++;

 // fill pull vector entry
 updates[cdc_index].V=fit_V;
}
else {
 updates[cdc_index].V=V;
}

// Flag that we used this hit
used_hits[cdc_index]=1;

// fill updates
updates[cdc_index].resi=res;
updates[cdc_index].d=d;
updates[cdc_index].delta=delta;
updates[cdc_index].S=S;
updates[cdc_index].C=C;
updates[cdc_index].tdrift=tdrift;
updates[cdc_index].ddrift=dmeas;
updates[cdc_index].s=29.98*trajectory[k].t; // assume beta=1
trajectory[k].id=cdc_index+1;
    }
    // move to next cdc hit
    if (cdc_index>0){
cdc_index--;

//New wire position
wire=cdchits[cdc_index]->wire;
if (VERBOSE>5) {
 jout << " Next Wire ring " << wire->ring << " straw " << wire->straw << " origin udir" << endl;
 wire->origin.Print(); wire->udir.Print();
}
origin=wire->origin;
z0=origin.z();
vz=wire->udir.z();
wdir=(1./vz)*wire->udir;
wirepos=origin+((trajectory[k].z-z0))*wdir;

// New doca^2
dx=S(state_x)-wirepos.x();
dy=S(state_y)-wirepos.y();
doca2=dx*dx+dy*dy;	

    }
    else more_hits=false;
  }
  
  old_doca2=doca2;
}
if (ndof<=4) return VALUE_OUT_OF_RANGE;

ndof-=4;

return NOERROR;
595}

597// Smooth the CDC only tracks
598jerror_t DTrackFitterStraightTrack::Smooth(vector<cdc_update_t>&cdc_updates){
unsigned int max=best_trajectory.size()-1;
DMatrix4x1 S=(best_trajectory[max].Skk);
DMatrix4x4 C=(best_trajectory[max].Ckk);
DMatrix4x4 JT=best_trajectory[max].J.Transpose();
DMatrix4x1 Ss=S;
DMatrix4x4 Cs=C;
DMatrix4x4 A,dC;
DMatrix1x4 H;  // Track projection matrix
DMatrix4x1 H_T; // Transpose of track projection matrix

const double d_EPS=1e-8;

for (unsigned int m=max-1;m>0;m--){
  if (best_trajectory[m].id>0){
    unsigned int id=best_trajectory[m].id-1;
    A=cdc_updates[id].C*JT*C.Invert();
    Ss=cdc_updates[id].S+A*(Ss-S);
    
    dC=A*(Cs-C)*A.Transpose();
    Cs=cdc_updates[id].C+dC;
    if (VERBOSE > 10) {
jout << " In Smoothing Step Using ID " << id << "/" << cdc_updates.size() << " for ring " << cdchits[id]->wire->ring << endl;
jout << " A cdc_updates[id].C Ss Cs " << endl;
A.Print(); cdc_updates[id].C.Print(); Ss.Print(); Cs.Print();
    }
    if(!Cs.IsPosDef()) {
if (VERBOSE) jout << "Cs is not PosDef!" << endl;
return VALUE_OUT_OF_RANGE;
    }

    const DCDCWire *wire=cdchits[id]->wire;
    DVector3 origin=wire->origin;
    double z0=origin.z();
    double vz=wire->udir.z();
    DVector3 wdir=(1./vz)*wire->udir;
    DVector3 wirepos=origin+(best_trajectory[m].z-z0)*wdir;
    // Position and direction from state vector
    double x=Ss(state_x);
    double y=Ss(state_y);
    double tx=Ss(state_tx);
    double ty=Ss(state_ty);
    
    DVector3 pos0(x,y,best_trajectory[m].z);
    DVector3 tdir(tx,ty,1.);
    
    // Find the true doca to the wire
    DVector3 diff=pos0-origin;
    double dx0=diff.x(),dy0=diff.y();
    double wdir_dot_diff=diff.Dot(wdir);
    double tdir_dot_diff=diff.Dot(tdir);
    double tdir_dot_wdir=tdir.Dot(wdir);
    double tdir2=tdir.Mag2();
    double wdir2=wdir.Mag2();
    double D=tdir2*wdir2-tdir_dot_wdir*tdir_dot_wdir;
    double N=tdir_dot_wdir*wdir_dot_diff-wdir2*tdir_dot_diff;
    double N1=tdir2*wdir_dot_diff-tdir_dot_wdir*tdir_dot_diff;
    double scale=1./D;
    double s=scale*N;
    double t=scale*N1;
    diff+=s*tdir-t*wdir;
    double d=diff.Mag()+d_EPS; // prevent division by zero
    double ddrift = cdc_updates[id].ddrift;
    
    double resi = ddrift - d;
    // Track projection
    
    {
double one_over_d=1./d;
double diffx=diff.x(),diffy=diff.y(),diffz=diff.z();
double wx=wdir.x(),wy=wdir.y();

double dN1dtx=2.*tx*wdir_dot_diff-wx*tdir_dot_diff-tdir_dot_wdir*dx0;
double dDdtx=2.*tx*wdir2-2.*tdir_dot_wdir*wx;
double dtdtx=scale*(dN1dtx-t*dDdtx);

double dN1dty=2.*ty*wdir_dot_diff-wy*tdir_dot_diff-tdir_dot_wdir*dy0;
double dDdty=2.*ty*wdir2-2.*tdir_dot_wdir*wy;
double dtdty=scale*(dN1dty-t*dDdty);

double dNdtx=wx*wdir_dot_diff-wdir2*dx0;
double dsdtx=scale*(dNdtx-s*dDdtx);

double dNdty=wy*wdir_dot_diff-wdir2*dy0;
double dsdty=scale*(dNdty-s*dDdty);

H(state_tx)=H_T(state_tx)
 =one_over_d*(diffx*(s+tx*dsdtx-wx*dtdtx)+diffy*(ty*dsdtx-wy*dtdtx)
       +diffz*(dsdtx-dtdtx));
H(state_ty)=H_T(state_ty)
 =one_over_d*(diffx*(tx*dsdty-wx*dtdty)+diffy*(s+ty*dsdty-wy*dtdty)
       +diffz*(dsdty-dtdty));

double dsdx=scale*(tdir_dot_wdir*wx-wdir2*tx);
double dtdx=scale*(tdir2*wx-tdir_dot_wdir*tx);
double dsdy=scale*(tdir_dot_wdir*wy-wdir2*ty);
double dtdy=scale*(tdir2*wy-tdir_dot_wdir*ty);

H(state_x)=H_T(state_x)
 =one_over_d*(diffx*(1.+dsdx*tx-dtdx*wx)+diffy*(dsdx*ty-dtdx*wy)
       +diffz*(dsdx-dtdx));
H(state_y)=H_T(state_y)
 =one_over_d*(diffx*(dsdy*tx-dtdy*wx)+diffy*(1.+dsdy*ty-dtdy*wy)
       +diffz*(dsdy-dtdy));
    }
    double V=cdc_updates[id].V;
  
    if (VERBOSE > 10) jout << " d " << d << " H*S " << H*S << endl;
    if (cdchits[id]->wire->ring==RING_TO_SKIP){
V=V+H*Cs*H_T;
    }
    else{
V=V-H*Cs*H_T;
    }      
    if (V<0) return VALUE_OUT_OF_RANGE;

    // Add the pull
    double myscale=1./sqrt(1.+tx*tx+ty*ty);
    double cosThetaRel=wire->udir.Dot(DVector3(myscale*tx,myscale*ty,myscale));
    pull_t thisPull(resi,sqrt(V),
      best_trajectory[m].t*SPEED_OF_LIGHT29.9792458,
      cdc_updates[id].tdrift,
      d,cdchits[id], NULL__null,
      diff.Phi(), //docaphi
      best_trajectory[m].z,cosThetaRel,
      cdc_updates[id].tdrift);
    
    // Derivatives for alignment
    double wtx=wire->udir.X(), wty=wire->udir.Y(), wtz=wire->udir.Z();
    double wx=wire->origin.X(), wy=wire->origin.Y(), wz=wire->origin.Z();
    
    double z=best_trajectory[m].z;
    double tx2=tx*tx, ty2=ty*ty;
    double wtx2=wtx*wtx, wty2=wty*wty, wtz2=wtz*wtz;
    double denom=(1 + ty2)*wtx2 + (1 + tx2)*wty2 - 2*ty*wty*wtz + (tx2 + ty2)*wtz2 - 2*tx*wtx*(ty*wty + wtz) +d_EPS;
    double denom2=denom*denom;
    double c1=-(wtx - tx*wtz)*(wy - y);
    double c2=wty*(wx - tx*wz - x + tx*z);
    double c3=ty*(-(wtz*wx) + wtx*wz + wtz*x - wtx*z);
    double dscale=0.5*(1./d);
    
    vector<double> derivatives(11);
    
    derivatives[CDCTrackD::dDOCAdOriginX]=dscale*(2*(wty - ty*wtz)*(c1 + c2 + c3))/denom;
    
    derivatives[CDCTrackD::dDOCAdOriginY]=dscale*(2*(-wtx + tx*wtz)*(c1 + c2 + c3))/denom;
    
    derivatives[CDCTrackD::dDOCAdOriginZ]=dscale*(2*(ty*wtx - tx*wty)*(c1 + c2 + c3))/denom;
    
    derivatives[CDCTrackD::dDOCAdDirX]=dscale*(2*(wty - ty*wtz)*(c1 + c2 + c3)*
				 (tx*(ty*wty + wtz)*(wx - x) + (wty - ty*wtz)*(-wy + y + ty*(wz - z)) +
				  wtx*(-((1 + ty2)*wx) + (1 + ty2)*x + tx*(ty*wy + wz - ty*y - z)) + tx2*(wty*(-wy + y) + wtz*(-wz + z))))/denom2;
    
    derivatives[CDCTrackD::dDOCAdDirY]=dscale*(-2*(wtx - tx*wtz)*(c1 + c2 + c3)*
				 (tx*(ty*wty + wtz)*(wx - x) + (wty - ty*wtz)*(-wy + y + ty*(wz - z)) +
				  wtx*(-((1 + ty2)*wx) + (1 + ty2)*x + tx*(ty*wy + wz - ty*y - z)) + tx2*(wty*(-wy + y) + wtz*(-wz + z))))/denom2;
    
    derivatives[CDCTrackD::dDOCAdDirZ]=dscale*(-2*(ty*wtx - tx*wty)*(c1 + c2 + c3)*
				 (-(tx*(ty*wty + wtz)*(wx - x)) + tx2*(wty*(wy - y) + wtz*(wz - z)) + (wty - ty*wtz)*(wy - y + ty*(-wz + z)) +
				  wtx*((1 + ty2)*wx - (1 + ty2)*x + tx*(-(ty*wy) - wz + ty*y + z))))/denom2;
    
    derivatives[CDCTrackD::dDOCAdS0]=-derivatives[CDCTrackD::dDOCAdOriginX];
    
    derivatives[CDCTrackD::dDOCAdS1]=-derivatives[CDCTrackD::dDOCAdOriginY];
    
    derivatives[CDCTrackD::dDOCAdS2]=dscale*(2*(wty - ty*wtz)*(-c1 - c2 - c3)*
			       (-(wtx*wtz*wx) - wty*wtz*wy + wtx2*wz + wty2*wz + wtx*wtz*x + wty*wtz*y - wtx2*z - wty2*z +
				tx*(wty2*(wx - x) + wtx*wty*(-wy + y) + wtz*(wtz*wx - wtx*wz - wtz*x + wtx*z)) +
				ty*(wtx*wty*(-wx + x) + wtx2*(wy - y) + wtz*(wtz*wy - wty*wz - wtz*y + wty*z))))/denom2;
    
    derivatives[CDCTrackD::dDOCAdS3]=dscale*(2*(wtx - tx*wtz)*(c1 + c2 + c3)*
			       (-(wtx*wtz*wx) - wty*wtz*wy + wtx2*wz + wty2*wz + wtx*wtz*x + wty*wtz*y - wtx2*z - wty2*z +
				tx*(wty2*(wx - x) + wtx*wty*(-wy + y) + wtz*(wtz*wx - wtx*wz - wtz*x + wtx*z)) +
				ty*(wtx*wty*(-wx + x) + wtx2*(wy - y) + wtz*(wtz*wy - wty*wz - wtz*y + wty*z))))/denom2;

    thisPull.AddTrackDerivatives(derivatives);
    pulls.push_back(thisPull);
  }
  else{
    A=best_trajectory[m].Ckk*JT*C.Invert();
    Ss=best_trajectory[m].Skk+A*(Ss-S);
    Cs=best_trajectory[m].Ckk+A*(Cs-C)*A.Transpose();
  }
  
  S=best_trajectory[m].Skk;
  C=best_trajectory[m].Ckk;
  JT=best_trajectory[m].J.Transpose();
}

return NOERROR;
788}

790//Reference trajectory for the track for cdc tracks
791jerror_t DTrackFitterStraightTrack::SetReferenceTrajectory(double t0,double z,
					   DMatrix4x1 &S,
					   const DCDCTrackHit *last_cdc,
					   double &dzsign){ 
DMatrix4x4 J(1.,0.,1.,0., 0.,1.,0.,1., 0.,0.,1.,0., 0.,0.,0.,1.);

double ds=0.5;
double t=t0;

// last y position of hit (approximate, using center of wire)
double last_y=last_cdc->wire->origin.y();
double last_r2=last_cdc->wire->origin.Perp2();
// Check that track points towards last wire, otherwise swap deltaz
DVector3 dir(S(state_tx),S(state_ty),dzsign);
if (fabs(dir.Theta() - M_PI_21.57079632679489661923) < 0.2) ds = 0.1;

//jout << "dPhi " << dphi << " theta " << dir.Theta() << endl;
double dz=dzsign*ds/sqrt(1.+S(state_tx)*S(state_tx)+S(state_ty)*S(state_ty));

if (VERBOSE) {
  if (COSMICS) jout << " Swimming Reference Trajectory last CDC y " << last_y << " dz "<< dz << endl;
  else jout << " Swimming Reference Trajectory last CDC r2 " << last_r2 << " dz "<< dz << endl;
  jout << " S" << endl; S.Print(); jout << "z= "<< z << endl;
  jout << "  Last CDC ring " << last_cdc->wire->ring << " straw " << last_cdc->wire->straw << endl;
}
unsigned int numsteps=0;
const unsigned int MAX_STEPS=5000;
bool done=false;
do{
  numsteps++;
  z+=dz;
  J(state_x,state_tx)=-dz;
  J(state_y,state_ty)=-dz;
  // Flight time: assume particle is moving at the speed of light
  t+=ds/29.98;
  //propagate the state to the next z position
  S(state_x)+=S(state_tx)*dz;
  S(state_y)+=S(state_ty)*dz;

  if (z > downstreamEndplate && dz < 0) continue;
  if (z < upstreamEndplate && dz > 0) continue;
  trajectory.push_front(trajectory_t(z,t,S,J,DMatrix4x1(),DMatrix4x4()));

  done=((z>downstreamEndplate) || (z<upstreamEndplate));
  if (COSMICS==false){
    double r2=S(state_x)*S(state_x)+S(state_y)*S(state_y);
    done |= (r2>last_r2);
  }
}while (!done && numsteps<MAX_STEPS);

if (VERBOSE)
  {
    if (VERBOSE > 10){
printf("Trajectory:\n");
for (unsigned int i=0;i<trajectory.size();i++){
 printf(" x %f y %f z %f\n",trajectory[i].S(state_x),
 trajectory[i].S(state_y),trajectory[i].z); 
}
    }
    else{
printf("%i step trajectory Begin/End:\n", numsteps);
printf(" x %f y %f z %f\n",trajectory[0].S(state_x), trajectory[0].S(state_y), trajectory[0].z);
if (trajectory.size() > 1) printf(" x %f y %f z %f\n",trajectory[trajectory.size()-1].S(state_x),
			  trajectory[trajectory.size()-1].S(state_y),trajectory[trajectory.size()-1].z);
    }
  }
if (trajectory.size()<2) return UNRECOVERABLE_ERROR;
return NOERROR;
859}

861// Routine for steering the fit of the track
862DTrackFitter::fit_status_t DTrackFitterStraightTrack::FitTrack(void){
DTrackFitter::fit_status_t status=DTrackFitter::kFitNotDone;
if (cdchits.size()+fdchits.size() < 4) return status;

// Initial guess for state vector
DVector3 input_pos=input_params.position();
DVector3 input_mom=input_params.momentum();
double pz=input_mom.z();
DMatrix4x1 S(input_pos.x(),input_pos.y(),input_mom.x()/pz,input_mom.y()/pz);

// Initial guess for covariance matrix
DMatrix4x4 C;
if (fit_type==kWireBased){
  C(state_x,state_x)=C(state_y,state_y)=4.0;     
  C(state_tx,state_tx)=C(state_ty,state_ty)=1.0;
}
else{
  C(state_x,state_x)=C(state_y,state_y)=1.;     
  C(state_tx,state_tx)=C(state_ty,state_ty)=0.01;
}

// Starting z-position and time
double z=input_pos.z();
double t0=input_params.t0();
// Variance in start time
switch(input_params.t0_detector()){
case SYS_TOF:
  mVarT0=0.01;
  break;
case SYS_CDC:
  mVarT0=25.;
  break;
case SYS_FDC:
  mVarT0=25.;
  break;
case SYS_BCAL:
  mVarT0=0.25;
  break;
case SYS_START:
  mVarT0=0.09;
  break;
default:
  mVarT0=0.;
  break;
}

// Chisq and ndof
chisq=1e16;
Ndof=0;

// Sort the CDC hits by radius or y (for cosmics)
if (cdchits.size()>0){
  if (COSMICS){
    stable_sort(cdchits.begin(),cdchits.end(),DTrackFitterStraightTrack_cdc_hit_cmp);
  }
  else{
    stable_sort(cdchits.begin(),cdchits.end(),DTrackFitterStraightTrack_cdc_hit_radius_cmp);
  }
}

if (fdchits.size()>0){
  // Sort FDC hits by z
  stable_sort(fdchits.begin(),fdchits.end(),DTrackFitterStraightTrack_fdc_hit_cmp);
  status=FitForwardTrack(t0,z,S,C,chisq,Ndof);
}
else if (cdchits.size()>0){
  double dzsign=(pz>0)?1.:-1.;
  status=FitCentralTrack(z,t0,dzsign,S,C,chisq,Ndof);
}

if (status==DTrackFitter::kFitSuccess){
  // Output fit results
  double tx=S(state_tx),ty=S(state_ty);
  double phi=atan2(ty,tx);
  double tanl=1./sqrt(tx*tx+ty*ty);
  // Check for tracks heading upstream
  if (cdchits_used_in_fit.size()>0){
    double phi_diff=phi-cdchits_used_in_fit[0]->wire->origin.Phi()-M_PI3.14159265358979323846;
    if (phi_diff<-M_PI3.14159265358979323846) phi_diff+=2.*M_PI3.14159265358979323846;
    if (phi_diff> M_PI3.14159265358979323846) phi_diff-=2.*M_PI3.14159265358979323846;
    if (fabs(phi_diff)<M_PI_21.57079632679489661923){
tanl*=-1;
phi+=M_PI3.14159265358979323846;
    }
  }
  double pt=5.0*cos(atan(tanl)); // arbitrary magnitude...
  DVector3 mom(pt*cos(phi),pt*sin(phi),pt*tanl);
  
  // Convert 4x4 covariance matrix to a TMatrixFSym for output
  TMatrixFSym errMatrix(5);
  for(unsigned int loc_i = 0; loc_i < 4; ++loc_i){
    for(unsigned int loc_j = 0; loc_j < 4; ++loc_j){
 errMatrix(loc_i, loc_j) = C(loc_i, loc_j);
    }
  }
  errMatrix(4,4)=1.;
  
  // Add 7x7 covariance matrix to output
  double sign=(mom.Theta()>M_PI_21.57079632679489661923)?-1.:1.;
  fit_params.setErrorMatrix(Get7x7ErrorMatrix(errMatrix,S,sign));

  DVector3 pos;
  if (COSMICS==false){
    DVector3 beamdir(0.,0.,1.);
    DVector3 beampos(0.,0.,65.);
    finder->FindDoca(z,S,beamdir,beampos,&pos);
    
    // Jacobian matrix
    double dz=pos.z()-z;
    DMatrix4x4 J(1.,0.,1.,0., 0.,1.,0.,1., 0.,0.,1.,0., 0.,0.,0.,1.);
    J(state_x,state_tx)=dz;
    J(state_y,state_ty)=dz;
    // Transform the matrix to the position of the doca
    C=J*C*J.Transpose();  
  }
  else{
    pos.SetXYZ(S(state_x),S(state_y),z);
  }

  fit_params.setForwardParmFlag(true);

  fit_params.setPosition(pos);
  fit_params.setMomentum(mom);

  // Add tracking covariance matrix to output
  auto locTrackingCovarianceMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();
  locTrackingCovarianceMatrix->ResizeTo(5, 5); 
  locTrackingCovarianceMatrix->Zero();
  for(unsigned int loc_i = 0; loc_i < 4; ++loc_i){
    for(unsigned int loc_j = 0; loc_j < 4; ++loc_j){
(*locTrackingCovarianceMatrix)(loc_i, loc_j) = C(loc_i, loc_j);
    }
  }
  (*locTrackingCovarianceMatrix)(4,4)=1.;
  fit_params.setTrackingErrorMatrix(locTrackingCovarianceMatrix);

  // Get extrapolations to other detectors
  GetExtrapolations(pos,mom);

}

return status;
1004}

1006// Locate a position in vector xx given x
1007unsigned int DTrackFitterStraightTrack::Locate(const vector<double>&xx,double x) const{
int n=xx.size();
if (x==xx[0]) return 0;
else if (x==xx[n-1]) return n-2;

int jl=-1;
int ju=n;
int ascnd=(xx[n-1]>=xx[0]);
while(ju-jl>1){
  int jm=(ju+jl)>>1;
  if ( (x>=xx[jm])==ascnd)
jl=jm;
  else
    ju=jm;
} 
return jl;
1023}

1025// Steering routine for fitting CDC-only tracks
1026DTrackFitter::fit_status_t DTrackFitterStraightTrack::FitCentralTrack(double &z0,double t0,
				double dzsign,
				DMatrix4x1 &Sbest,
				DMatrix4x4 &Cbest,
				double &chi2_best,
				int &ndof_best){
// State vector and covariance matrix
DMatrix4x1 S(Sbest);
DMatrix4x4 C(Cbest),C0(Cbest);

double chi2=chi2_best;
int ndof=ndof_best;

unsigned int numhits=cdchits.size();
unsigned int maxindex=numhits-1;

// vectors of indexes to cdc hits used in the fit
vector<int> used_cdc_hits(numhits);
vector<int> used_cdc_hits_best_fit(numhits);

// vectors of residual information 
vector<cdc_update_t>updates(numhits);
vector<cdc_update_t>best_updates(numhits);

// Rest deque for "best" trajectory
best_trajectory.clear();

// Perform the fit
int iter=0;
for(iter=0;iter<20;iter++){
  if (VERBOSE) jout << " Performing Pass iter " << iter << endl;
     
  trajectory.clear();
  if (SetReferenceTrajectory(t0,z0,S,cdchits[maxindex],dzsign)
!=NOERROR) break;
  
  if (VERBOSE) jout << " Reference Trajectory Set " << endl;
  C=C0;
  if (KalmanFilter(S,C,used_cdc_hits,updates,chi2,ndof,iter)!=NOERROR) break;
  if (VERBOSE) jout << " Filter returns NOERROR" << endl;
  if (iter>0 && (fabs(chi2_best-chi2)<0.1 || chi2>chi2_best)) break;  
  
  // Save the current state and covariance matrixes
  Cbest=C;
  Sbest=S;

  // Save the current used hit and trajectory information
  best_trajectory.assign(trajectory.begin(),trajectory.end());
  used_cdc_hits_best_fit.assign(used_cdc_hits.begin(),used_cdc_hits.end());
  best_updates.assign(updates.begin(),updates.end());

  chi2_best=chi2; 
  ndof_best=ndof;
}
if (iter==0) return DTrackFitter::kFitFailed;

//Final z position (closest to beam line)
z0=best_trajectory[best_trajectory.size()-1].z;

//Run the smoother 
if (Smooth(best_updates) == NOERROR) IsSmoothed=true;

// output list of cdc hits used in the fit
cdchits_used_in_fit.clear();
for (unsigned int m=0;m<used_cdc_hits_best_fit.size();m++){
  if (used_cdc_hits_best_fit[m]){
    cdchits_used_in_fit.push_back(cdchits[m]);
  }
}

return DTrackFitter::kFitSuccess;
1097}


1100//----------------------------------------------------
1101// FDC fitting routines
1102//----------------------------------------------------

1104// parametrization of time-to-distance for FDC
1105double DTrackFitterStraightTrack::fdc_drift_distance(double time) const {
if (time<0.) return 0.;
double d=0.; 
double tsq=time*time;
double t_high=DRIFT_FUNC_PARMS[4];

if (time<t_high){
  d=DRIFT_FUNC_PARMS[0]*sqrt(time)+DRIFT_FUNC_PARMS[1]*time
    +DRIFT_FUNC_PARMS[2]*tsq+DRIFT_FUNC_PARMS[3]*tsq*time;
}
else{
  double t_high_sq=t_high*t_high;
  d=DRIFT_FUNC_PARMS[0]*sqrt(t_high)+DRIFT_FUNC_PARMS[1]*t_high
    +DRIFT_FUNC_PARMS[2]*t_high_sq+DRIFT_FUNC_PARMS[3]*t_high_sq*t_high;
  d+=DRIFT_FUNC_PARMS[5]*(time-t_high);
}
  
return d;

1124}

1126// Crude approximation for the variance in drift distance due to smearing
1127double DTrackFitterStraightTrack::fdc_drift_variance(double t) const{
 //return FDC_ANODE_VARIANCE;
 if (t<5.) t=5.;
 double sigma=DRIFT_RES_PARMS[0]/(t+1.)+DRIFT_RES_PARMS[1]+DRIFT_RES_PARMS[2]*t*t;

 return sigma*sigma;
1133}

1135// Steering routine for the kalman filter
1136DTrackFitter::fit_status_t
1137DTrackFitterStraightTrack::FitForwardTrack(double t0,double &start_z,
    DMatrix4x1 &Sbest,DMatrix4x4 &Cbest,double &chi2_best,int &ndof_best){
// State vector and covariance matrix
DMatrix4x1 S(Sbest);
DMatrix4x4 C(Cbest),C0(Cbest);

// vectors of indexes to fdc hits used in the fit
unsigned int numfdchits=fdchits.size();
vector<int> used_fdc_hits(numfdchits);
vector<int> used_fdc_hits_best_fit(numfdchits);
// vectors of indexes to cdc hits used in the fit
unsigned int numcdchits=cdchits.size();
vector<int> used_cdc_hits(numcdchits);
vector<int> used_cdc_hits_best_fit(numcdchits);

// vectors of residual information 
vector<fdc_update_t>updates(numfdchits);
vector<fdc_update_t>best_updates(numfdchits);
vector<cdc_update_t>cdc_updates(numcdchits);
vector<cdc_update_t>best_cdc_updates(numcdchits);

vector<const DCDCTrackHit *> matchedCDCHits;

// Chi-squared and degrees of freedom
double chi2=chi2_best;
int ndof=ndof_best;

// Rest deque for "best" trajectory
best_trajectory.clear();

unsigned iter=0;
// First pass
for(iter=0;iter<20;iter++){
  trajectory.clear();
  if (SetReferenceTrajectory(t0,start_z,S)!=NOERROR) break;
  
  C=C0;
  if (KalmanFilter(S,C,used_fdc_hits,used_cdc_hits,updates,cdc_updates,
     chi2,ndof)!=NOERROR) break;
  
  // printf(" == iter %d =====chi2 %f ndof %d \n",iter,chi2,ndof);
  if (iter>0 && (chi2>chi2_best || fabs(chi2_best-chi2)<0.1)) break;  
  
  // Save the current state and covariance matrixes
  Cbest=C;
  Sbest=S;
  
  // Save the current used hit and trajectory information
  best_trajectory.assign(trajectory.begin(),trajectory.end());
  used_cdc_hits_best_fit.assign(used_cdc_hits.begin(),used_cdc_hits.end());
  used_fdc_hits_best_fit.assign(used_fdc_hits.begin(),used_fdc_hits.end());
  best_updates.assign(updates.begin(),updates.end());
  best_cdc_updates.assign(cdc_updates.begin(),cdc_updates.end());
  
  chi2_best=chi2; 
  ndof_best=ndof;
}
if (iter==0) return DTrackFitter::kFitFailed;

//Final z position (closest to beam line)
start_z=best_trajectory[best_trajectory.size()-1].z;

//Run the smoother
if (Smooth(best_updates,best_cdc_updates) == NOERROR) IsSmoothed=true;

// output list of hits used in the fit
cdchits_used_in_fit.clear();
for (unsigned int m=0;m<used_cdc_hits_best_fit.size();m++){
  if (used_cdc_hits_best_fit[m]){
    cdchits_used_in_fit.push_back(cdchits[m]);
  }
}  
fdchits_used_in_fit.clear();
for (unsigned int m=0;m<used_fdc_hits_best_fit.size();m++){
  if (used_fdc_hits_best_fit[m]){
    fdchits_used_in_fit.push_back(fdchits[m]);
  }
}


return DTrackFitter::kFitSuccess;
1218}


1221// Reference trajectory for the track
1222jerror_t 
1223DTrackFitterStraightTrack::SetReferenceTrajectory(double t0,double z,
				  DMatrix4x1 &S){
const double EPS=1e-3;

// Jacobian matrix 
DMatrix4x4 J(1.,0.,1.,0., 0.,1.,0.,1., 0.,0.,1.,0., 0.,0.,0.,1.);

double dz=1.1;
double t=t0;
trajectory.push_front(trajectory_t(z,t,S,J,DMatrix4x1(),DMatrix4x4()));

double zhit=z;
double old_zhit=z;
for (unsigned int i=0;i<fdchits.size();i++){  
  zhit=fdchits[i]->wire->origin.z();
  dz=1.1;
  
  if (fabs(zhit-old_zhit)<EPS && i>0){
    trajectory[0].numhits++;
    continue;
  }
  // propagate until we would step beyond the FDC hit plane
  bool done=false;
  while (!done){	    
    double new_z=z+dz;	      
    
    if (new_z>zhit){
dz=zhit-z;
new_z=zhit;
done=true;
    }
    J(state_x,state_tx)=-dz;
    J(state_y,state_ty)=-dz;
    // Flight time: assume particle is moving at the speed of light
    t+=dz*sqrt(1+S(state_tx)*S(state_tx)+S(state_ty)*S(state_ty))/29.98;
    //propagate the state to the next z position
    S(state_x)+=S(state_tx)*dz;
    S(state_y)+=S(state_ty)*dz;
    
    trajectory.push_front(trajectory_t(new_z,t,S,J,DMatrix4x1(),
			 DMatrix4x4())); 
    if (done){
trajectory[0].id=i+1;
trajectory[0].numhits=1;
    }
    
    z=new_z;
  }	   
  old_zhit=zhit;
}
// One last step
dz=1.1;
J(state_x,state_tx)=-dz;
J(state_y,state_ty)=-dz;

 // Flight time: assume particle is moving at the speed of light
t+=dz*sqrt(1+S(state_tx)*S(state_tx)+S(state_ty)*S(state_ty))/29.98;

//propagate the state to the next z position
S(state_x)+=S(state_tx)*dz;
S(state_y)+=S(state_ty)*dz;
trajectory.push_front(trajectory_t(z+dz,t,S,J,DMatrix4x1(),DMatrix4x4()));

if (false)
 {
   printf("Trajectory:\n");
   for (unsigned int i=0;i<trajectory.size();i++){
     printf(" x %f y %f z %f first hit %d num in layer %d\n",trajectory[i].S(state_x),
     trajectory[i].S(state_y),trajectory[i].z,trajectory[i].id,
     trajectory[i].numhits); 
   }
 }

return NOERROR;
1297}

1299// Perform Kalman Filter for the current trajectory
1300jerror_t DTrackFitterStraightTrack::KalmanFilter(DMatrix4x1 &S,DMatrix4x4 &C,
				 vector<int>&used_fdc_hits, 
				 vector<int>&used_cdc_hits,
				 vector<fdc_update_t>&updates,
				 vector<cdc_update_t>&cdc_updates,
				 double &chi2,int &ndof){
DMatrix2x4 H;  // Track projection matrix
DMatrix4x2 H_T; // Transpose of track projection matrix 
DMatrix4x2 K;  // Kalman gain matrix
double VarMs=0.001; // kludge for material
DMatrix2x2 V(0.0833,0.,0.,0.000256+VarMs);  // Measurement variance 
DMatrix2x2 Vtemp,InvV;
DMatrix2x1 Mdiff;
DMatrix4x4 I; // identity matrix
DMatrix4x4 J; // Jacobian matrix
DMatrix4x1 S0; // State vector from reference trajectory

DMatrix1x4 H_CDC;  // Track projection matrix
DMatrix4x1 H_T_CDC; // Transpose of track projection matrix
DMatrix4x1 K_CDC;  // Kalman gain matrix
double V_CDC;

const double d_EPS=1e-8;

// Zero out the vectors of used hit flags
for (unsigned int i=0;i<used_fdc_hits.size();i++) used_fdc_hits[i]=0; 
for (unsigned int i=0;i<used_cdc_hits.size();i++) used_cdc_hits[i]=0;

//Initialize chi2 and ndof
chi2=0.;
ndof=0;

// Save the starting values for C and S in the deque
trajectory[0].Skk=S;
trajectory[0].Ckk=C;

// Loop over all steps in the trajectory
S0=trajectory[0].S;
J=trajectory[0].J;

// CDC index and wire position variables
bool more_hits = cdchits.size() == 0 ? false: true;
bool firstCDCStep=true;
unsigned int cdc_index=0;
const DCDCWire *wire;
DVector3 origin,wdir,wirepos;
double doca2=0.0, old_doca2=0.0;
if(more_hits){
  cdc_index=cdchits.size()-1;
  wire=cdchits[cdc_index]->wire;
  origin=wire->origin;
  double vz=wire->udir.z();
  if (VERBOSE) jout << " Additional CDC Hits in FDC track Starting in Ring " << wire->ring << endl;
  wdir=(1./vz)*wire->udir;
}

for (unsigned int k=1;k<trajectory.size();k++){
  if (C(0,0)<=0. || C(1,1)<=0. || C(2,2)<=0. || C(3,3)<=0.)
    return UNRECOVERABLE_ERROR;
  
  // Propagate the state and covariance matrix forward in z
  S=trajectory[k].S+J*(S-S0);
  C=J*C*J.Transpose();
  
  // Save the current state and covariance matrix in the deque
  trajectory[k].Skk=S;
  trajectory[k].Ckk=C;
   
  // Save S and J for the next step
  S0=trajectory[k].S;
  J=trajectory[k].J;
  
  // Correct S and C for the hit 
  if (trajectory[k].id>0){
    unsigned int id=trajectory[k].id-1;
    
    double cospsi=cos(fdchits[id]->wire->angle);
    double sinpsi=sin(fdchits[id]->wire->angle);

    // Small angle alignment correction
    double x = S(state_x) + fdchits[id]->wire->angles.Z()*S(state_y);
    double y = S(state_y) - fdchits[id]->wire->angles.Z()*S(state_x);
    //tz = 1. + my_fdchits[id]->phiY*tx - my_fdchits[id]->phiX*ty;
    double tx = (S(state_tx) + fdchits[id]->wire->angles.Z()*S(state_ty) - fdchits[id]->wire->angles.Y());
    double ty = (S(state_ty) - fdchits[id]->wire->angles.Z()*S(state_tx) + fdchits[id]->wire->angles.X());
    
    if (std::isnan(x) || std::isnan(y)) return UNRECOVERABLE_ERROR;
    
    // x,y and tx,ty in local coordinate system	
    // To transform from (x,y) to (u,v), need to do a rotation:
    //   u = x*cos(psi)-y*sin(psi)
    //   v = y*cos(psi)+x*sin(psi)
    // (without alignment offsets)
    double vpred_wire_plane=y*cospsi+x*sinpsi;
    double upred_wire_plane=x*cospsi-y*sinpsi;
    double tu=tx*cospsi-ty*sinpsi;
    double tv=tx*sinpsi+ty*cospsi;
    
    // Variables for angle of incidence with respect to the z-direction in
    // the u-z plane
    double alpha=atan(tu);
    double cosalpha=cos(alpha);
    double cos2_alpha=cosalpha*cosalpha;
    double sinalpha=sin(alpha);
    double sin2_alpha=sinalpha*sinalpha;
    double cos2_alpha_minus_sin2_alpha=cos2_alpha-sin2_alpha;
    
    // Difference between measurement and projection
    for (int m=trajectory[k].numhits-1;m>=0;m--){
unsigned int my_id=id+m;
double uwire=fdchits[my_id]->w;
// (signed) distance of closest approach to wire
double du=upred_wire_plane-uwire;
double doca=du*cosalpha;

// Predicted avalanche position along the wire
double vpred=vpred_wire_plane;

// Measured position of hit along wire
double v=fdchits[my_id]->s; 

// Difference between measurements and predictions
double drift=0.; // assume hit at wire position
if (fit_type==kTimeBased){
 double drift_time=fdchits[my_id]->time-trajectory[k].t; 
 drift=(du>0.0?1.:-1.)*fdc_drift_distance(drift_time);	 
 V(0,0)=fdc_drift_variance(drift_time)+VarMs;
}
Mdiff(0)=drift-doca;
Mdiff(1)=v-vpred;

// Matrix for transforming from state-vector space to measurement space
H_T(state_x,0)=cospsi*cosalpha;
H_T(state_y,0)=-sinpsi*cosalpha;
double temp=-du*sinalpha*cos2_alpha;
H_T(state_tx,0)=cospsi*temp;
H_T(state_ty,0)=-sinpsi*temp;
double temp2=cosalpha*sinalpha*tv;
H_T(state_x,1)=sinpsi-temp2*cospsi;
H_T(state_y,1)=cospsi+temp2*sinpsi;
double temp4=sinalpha*doca;
double temp5=tv*cos2_alpha*du*cos2_alpha_minus_sin2_alpha;
H_T(state_tx,1)=-sinpsi*temp4-cospsi*temp5;
H_T(state_ty,1)=-cospsi*temp4+sinpsi*temp5;

// Matrix transpose H_T -> H
H(0,state_x)=H_T(state_x,0);
H(0,state_y)=H_T(state_y,0);
H(0,state_tx)=H_T(state_tx,0);
H(0,state_ty)=H_T(state_ty,0);
H(1,state_x)=H_T(state_x,1);
H(1,state_y)=H_T(state_y,1);
H(1,state_tx)=H_T(state_tx,1);
H(1,state_ty)=H_T(state_ty,1);

// Variance for this hit
InvV=(V+H*C*H_T).Invert();

// Compute Kalman gain matrix
K=(C*H_T)*InvV;

if (fdchits[my_id]->wire->layer!=PLANE_TO_SKIP){
 /*
   if(DEBUG_HISTS){
   hFDCOccTrkFit->Fill(fdchits[my_id]->wire->layer);
   }
 */
 // Update the state vector 
 S+=K*Mdiff;
 if(VERBOSE) S.Print();
 // Update state vector covariance matrix
 C=C-K*(H*C);    
 
 // Update the filtered measurement covariane matrix and put results in 
 // update vector
 DMatrix2x2 RC=V-H*C*H_T;
 DMatrix2x1 res=Mdiff-H*K*Mdiff;
 
 chi2+=RC.Chi2(res);
 ndof+=2;
  
 // fill pull vector entries
 updates[my_id].V=RC;
}
else{
 updates[my_id].V=V;
}

used_fdc_hits[my_id]=1;

// fill pull vector
updates[my_id].d=doca;
updates[my_id].S=S;
updates[my_id].C=C;
updates[my_id].tdrift=fdchits[my_id]->time-trajectory[k].t;
updates[my_id].s=29.98*trajectory[k].t; // assume beta=1
    } 
  }
  
  if (more_hits && trajectory[k].z < downstreamEndplate){
    // Position along wire
    double z0=origin.Z();
    wirepos=origin+(trajectory[k].z-z0)*wdir;
    
    // New doca^2
    double dx=S(state_x)-wirepos.X();
    double dy=S(state_y)-wirepos.Y();
    doca2=dx*dx+dy*dy;
    if (VERBOSE > 10) jout<< "At Position " << S(state_x) << " " << S(state_y) << " " << trajectory[k].z << " doca2 " << doca2 << endl;
    
    if (doca2>old_doca2 && more_hits && !firstCDCStep){

// zero-position and direction of line describing particle trajectory
double tx=S(state_tx),ty=S(state_ty);
DVector3 pos0(S(state_x),S(state_y),trajectory[k].z);
DVector3 tdir(tx,ty,1.);

// Find the true doca to the wire
DVector3 diff=pos0-origin;
double dx0=diff.x(),dy0=diff.y();
double wdir_dot_diff=diff.Dot(wdir);
double tdir_dot_diff=diff.Dot(tdir);
double tdir_dot_wdir=tdir.Dot(wdir);
double tdir2=tdir.Mag2();
double wdir2=wdir.Mag2();
double D=tdir2*wdir2-tdir_dot_wdir*tdir_dot_wdir;
double N=tdir_dot_wdir*wdir_dot_diff-wdir2*tdir_dot_diff;
double N1=tdir2*wdir_dot_diff-tdir_dot_wdir*tdir_dot_diff;
double scale=1./D;
double s=scale*N;
double t=scale*N1;
diff+=s*tdir-t*wdir;
double d=diff.Mag()+d_EPS; // prevent division by zero

// The next measurement	
double dmeas=0.39,delta=0.;  
double tdrift=cdchits[cdc_index]->tdrift-trajectory[k].t;	 
if (fit_type==kTimeBased){	 
 double phi_d=diff.Phi();
 double dphi=phi_d-origin.Phi();
 while (dphi>M_PI3.14159265358979323846) dphi-=2*M_PI3.14159265358979323846;
 while (dphi<-M_PI3.14159265358979323846) dphi+=2*M_PI3.14159265358979323846;
 
 int ring_index=cdchits[cdc_index]->wire->ring-1;
 int straw_index=cdchits[cdc_index]->wire->straw-1;
 double dz=t*wdir.z();
 double delta=max_sag[ring_index][straw_index]*(1.-dz*dz/5625.)
   *cos(phi_d+sag_phi_offset[ring_index][straw_index]);
 CDCDriftParameters(dphi,delta,tdrift,dmeas,V_CDC);
}

// residual
double res=dmeas-d;
if (VERBOSE>5) jout << " Residual " << res << endl;
// Track projection
double one_over_d=1./d;
double diffx=diff.x(),diffy=diff.y(),diffz=diff.z();
double wx=wdir.x(),wy=wdir.y();

double dN1dtx=2.*tx*wdir_dot_diff-wx*tdir_dot_diff-tdir_dot_wdir*dx0;
double dDdtx=2.*tx*wdir2-2.*tdir_dot_wdir*wx;
double dtdtx=scale*(dN1dtx-t*dDdtx);

double dN1dty=2.*ty*wdir_dot_diff-wy*tdir_dot_diff-tdir_dot_wdir*dy0;
double dDdty=2.*ty*wdir2-2.*tdir_dot_wdir*wy;
double dtdty=scale*(dN1dty-t*dDdty);

double dNdtx=wx*wdir_dot_diff-wdir2*dx0;
double dsdtx=scale*(dNdtx-s*dDdtx);

double dNdty=wy*wdir_dot_diff-wdir2*dy0;
double dsdty=scale*(dNdty-s*dDdty);

H_CDC(state_tx)=H_T_CDC(state_tx)
 =one_over_d*(diffx*(s+tx*dsdtx-wx*dtdtx)+diffy*(ty*dsdtx-wy*dtdtx)
       +diffz*(dsdtx-dtdtx));
H_CDC(state_ty)=H_T_CDC(state_ty)
 =one_over_d*(diffx*(tx*dsdty-wx*dtdty)+diffy*(s+ty*dsdty-wy*dtdty)
       +diffz*(dsdty-dtdty));

double dsdx=scale*(tdir_dot_wdir*wx-wdir2*tx);
double dtdx=scale*(tdir2*wx-tdir_dot_wdir*tx);
double dsdy=scale*(tdir_dot_wdir*wy-wdir2*ty);
double dtdy=scale*(tdir2*wy-tdir_dot_wdir*ty);

H_CDC(state_x)=H_T_CDC(state_x)
 =one_over_d*(diffx*(1.+dsdx*tx-dtdx*wx)+diffy*(dsdx*ty-dtdx*wy)
       +diffz*(dsdx-dtdx));
H_CDC(state_y)=H_T_CDC(state_y)
 =one_over_d*(diffx*(dsdy*tx-dtdy*wx)+diffy*(1.+dsdy*ty-dtdy*wy)
       +diffz*(dsdy-dtdy));

double InvV=1./(V_CDC+H_CDC*C*H_T_CDC);

// Check how far this hit is from the projection
double chi2check=res*res*InvV;
if (chi2check < CHI2CUT || DO_PRUNING == 0){
 if (VERBOSE) jout << "CDC Hit Added to FDC track " << endl;
 if (cdchits[cdc_index]->wire->ring!=RING_TO_SKIP){

   // Compute Kalman gain matrix
   K_CDC=InvV*(C*H_T_CDC);
   // Update state vector covariance matrix
   DMatrix4x4 Ctest=C-K_CDC*(H_CDC*C);
   
   //C.Print();
   //K.Print();
   //Ctest.Print();
   
   // Check that Ctest is positive definite
   if (!Ctest.IsPosDef()) return VALUE_OUT_OF_RANGE;
   C=Ctest;
   if(VERBOSE>10) C.Print();
   // Update the state vector
   //S=S+res*K;
   S+=res*K_CDC;
   if(VERBOSE) {jout << "traj[z]=" << trajectory[k].z<< endl; S.Print();} 
   
   // Compute new residual
   //d=finder->FindDoca(trajectory[k].z,S,wdir,origin);
   res=res-H_CDC*K_CDC*res;
   
   // Update chi2
   double fit_V=V_CDC-H_CDC*C*H_T_CDC;
   chi2+=res*res/fit_V;
   ndof++;
 
   // fill pull vector entry
   cdc_updates[cdc_index].V=fit_V;
 }
 else {
   cdc_updates[cdc_index].V=V_CDC;
 }

 // fill updates
 cdc_updates[cdc_index].resi=res;
 cdc_updates[cdc_index].d=d;
 cdc_updates[cdc_index].delta=delta;
 cdc_updates[cdc_index].S=S;
 cdc_updates[cdc_index].C=C;
 cdc_updates[cdc_index].tdrift=tdrift;
 cdc_updates[cdc_index].ddrift=dmeas;
 cdc_updates[cdc_index].s=29.98*trajectory[k].t; // assume beta=1
 trajectory[k].id=cdc_index+1000;
 
 used_cdc_hits[cdc_index]=1;
}
// move to next cdc hit
if (cdc_index>0){
 cdc_index--;
 
 //New wire position
 wire=cdchits[cdc_index]->wire;
 if (VERBOSE>5) jout << " Next Wire ring " << wire->ring << endl;
 origin=wire->origin;
 double vz=wire->udir.z();
 wdir=(1./vz)*wire->udir;
 wirepos=origin+((trajectory[k].z-z0))*wdir;
 
 // New doca^2
 dx=S(state_x)-wirepos.x();
 dy=S(state_y)-wirepos.y();
 doca2=dx*dx+dy*dy;
 
}
else more_hits=false;
    }
    firstCDCStep=false;
    old_doca2=doca2;
  }
}

ndof-=4;

return NOERROR;
1675}


1678// Smoothing algorithm for the forward trajectory.  Updates the state vector
1679// at each step (going in the reverse direction to the filter) based on the 
1680// information from all the steps and outputs the state vector at the
1681// outermost step.

1683jerror_t 
1684DTrackFitterStraightTrack::Smooth(vector<fdc_update_t>&fdc_updates,
    vector<cdc_update_t>&cdc_updates){ 
unsigned int max=best_trajectory.size()-1;
DMatrix4x1 S=(best_trajectory[max].Skk);
DMatrix4x4 C=(best_trajectory[max].Ckk);
DMatrix4x4 JT=best_trajectory[max].J.Transpose();
DMatrix4x1 Ss=S;
DMatrix4x4 Cs=C;
DMatrix4x4 A,dC;

const double d_EPS=1e-8;

for (unsigned int m=max-1;m>0;m--){
  if (best_trajectory[m].id>0 && best_trajectory[m].id<1000){ // FDC Hit
    unsigned int id=best_trajectory[m].id-1;
    A=fdc_updates[id].C*JT*C.Invert();
    Ss=fdc_updates[id].S+A*(Ss-S);

    dC=A*(Cs-C)*A.Transpose();
    Cs=fdc_updates[id].C+dC;

    double cosa=cos(fdchits[id]->wire->angle);
    double cos2a=cos(2*fdchits[id]->wire->angle);
    double sina=sin(fdchits[id]->wire->angle);
    double u=fdchits[id]->w;
    double v=fdchits[id]->s;

    // Small angle alignment correction
    double x = Ss(state_x) + fdchits[id]->wire->angles.Z() * Ss(state_y);
    double y = Ss(state_y) - fdchits[id]->wire->angles.Z() * Ss(state_x);
    // tz = 1.0 + my_fdchits[id]->phiY * tx - my_fdchits[id]->phiX * ty;
    double tx = Ss(state_tx) + fdchits[id]->wire->angles.Z() * Ss(state_ty) - fdchits[id]->wire->angles.Y();
    double ty = Ss(state_ty) - fdchits[id]->wire->angles.Z() * Ss(state_tx) + fdchits[id]->wire->angles.X();

    // Projected position along the wire 
    double vpred=x*sina+y*cosa;
    
    // Projected position in the plane of the wires transverse to the wires
    double upred=x*cosa-y*sina;

    // Direction tangent in the u-z plane
    double tu=tx*cosa-ty*sina;
    double alpha=atan(tu);
    double cosalpha=cos(alpha);
    //double cosalpha2=cosalpha*cosalpha;
    double sinalpha=sin(alpha);

    // (signed) distance of closest approach to wire
    double du=upred-u;
    double doca=du*cosalpha;
    // Difference between measurement and projection for the cathodes
    double tv=tx*sina+ty*cosa;
    double resi_c=v-vpred;

    // Difference between measurement and projection perpendicular to the wire
    double drift = 0.0;  // assume hit at wire position
    int left_right = -999;
    double drift_time = fdc_updates[id].tdrift;
    if (fit_type == kTimeBased) {
      drift = (du > 0.0 ? 1.0 : -1.0) * fdc_drift_distance(drift_time);
      left_right = (du > 0.0 ? +1 : -1);
    }
    double resi_a = drift - doca;

    // Variance from filter step
    DMatrix2x2 V=fdc_updates[id].V;
    // Compute projection matrix and find the variance for the residual
    DMatrix4x2 H_T;
    double temp2=-tv*sinalpha;
    H_T(state_x,1)=sina+cosa*cosalpha*temp2;	
    H_T(state_y,1)=cosa-sina*cosalpha*temp2;	
    
    double cos2_minus_sin2=cosalpha*cosalpha-sinalpha*sinalpha;
    double doca_cosalpha=doca*cosalpha;
    H_T(state_tx,1)=-doca_cosalpha*(tu*sina+tv*cosa*cos2_minus_sin2);
    H_T(state_ty,1)=-doca_cosalpha*(tu*cosa-tv*sina*cos2_minus_sin2);
    
    H_T(state_x,0)=cosa*cosalpha;
    H_T(state_y,0)=-sina*cosalpha;
    double one_plus_tu2=1.+tu*tu;
    double factor=du*tu/sqrt(one_plus_tu2)/one_plus_tu2;
    H_T(state_ty,0)=sina*factor;
    H_T(state_tx,0)=-cosa*factor;

    // Matrix transpose H_T -> H
    DMatrix2x4 H;
    H(0,state_x)=H_T(state_x,0);
    H(0,state_y)=H_T(state_y,0);
    H(0,state_tx)=H_T(state_tx,0);
    H(0,state_ty)=H_T(state_ty,0);
    H(1,state_x)=H_T(state_x,1);
    H(1,state_y)=H_T(state_y,1);
    H(1,state_tx)=H_T(state_tx,1);
    H(1,state_ty)=H_T(state_ty,1);

    if (fdchits[id]->wire->layer == PLANE_TO_SKIP) {
      // V += Cs.SandwichMultiply(H_T);
      V = V + H * Cs * H_T;
    } else {
      // V -= dC.SandwichMultiply(H_T);

      // R. Fruehwirth, Nucl. Instrum. Methods Phys. Res. A 262, 444 (1987).
      // Eq. (9)
      // The following V (lhs) corresponds to R^n_k in the paper.
      // dC corresponds to 'A_k * (C^n_{k+1} - C^k_{k+1}) * A_k^T' in the paper.
      V = V - H * dC * H_T;
    }

    /*
    if(DEBUG_HISTS){
hFDCOccTrkSmooth->Fill(fdchits[id]->wire->layer);
    }
    */
    // Implement derivatives wrt track parameters needed for millepede alignment
    // Add the pull
    double scale=1./sqrt(1.+tx*tx+ty*ty);
    double cosThetaRel=fdchits[id]->wire->udir.Dot(DVector3(scale*tx,scale*ty,scale));
    DTrackFitter::pull_t thisPull(resi_a,sqrt(V(0,0)),
		    best_trajectory[m].t*SPEED_OF_LIGHT29.9792458,
		    fdc_updates[id].tdrift,
		    fdc_updates[id].d,
		    NULL__null,fdchits[id],
		    0.0, //docaphi
		    best_trajectory[m].z,cosThetaRel, 
		    0.0, //tcorr
		    resi_c, sqrt(V(1,1))
		    );
    thisPull.left_right = left_right;

    if (fdchits[id]->wire->layer!=PLANE_TO_SKIP){
vector<double> derivatives;
derivatives.resize(FDCTrackD::size);

//dDOCAW/dDeltaX
derivatives[FDCTrackD::dDOCAW_dDeltaX] = -(1/sqrt(1 + pow(tx*cosa - ty*sina,2)));

//dDOCAW/dDeltaZ
derivatives[FDCTrackD::dDOCAW_dDeltaZ] = (tx*cosa - ty*sina)/sqrt(1 + pow(tx*cosa - ty*sina,2));

//dDOCAW/ddeltaPhiX
derivatives[FDCTrackD::dDOCAW_dDeltaPhiX] = (sina*(-(tx*cosa) + ty*sina)*(u - x*cosa + y*sina))/pow(1 + pow(tx*cosa - ty*sina,2),1.5);

//dDOCAW/ddeltaphiY
derivatives[FDCTrackD::dDOCAW_dDeltaPhiY] = (cosa*(tx*cosa - ty*sina)*(-u + x*cosa - y*sina))/pow(1 + pow(tx*cosa - ty*sina,2),1.5);

//dDOCAW/ddeltaphiZ
derivatives[FDCTrackD::dDOCAW_dDeltaPhiZ] = (tx*ty*u*cos2a + (x + pow(ty,2)*x - tx*ty*y)*sina + 
				     cosa*(-(tx*ty*x) + y + pow(tx,2)*y + (pow(tx,2) - pow(ty,2))*u*sina))/
 pow(1 + pow(tx*cosa - ty*sina,2),1.5);

// dDOCAW/dx
derivatives[FDCTrackD::dDOCAW_dx] = cosa/sqrt(1 + pow(tx*cosa - ty*sina,2));

// dDOCAW/dy
derivatives[FDCTrackD::dDOCAW_dy] = -(sina/sqrt(1 + pow(tx*cosa - ty*sina,2)));

// dDOCAW/dtx
derivatives[FDCTrackD::dDOCAW_dtx] = -((cosa*(tx*cosa - ty*sina)*(-u + x*cosa - y*sina))/pow(1 + pow(tx*cosa - ty*sina,2),1.5));

// dDOCAW/dty
derivatives[FDCTrackD::dDOCAW_dty] = (sina*(-(tx*cosa) + ty*sina)*(u - x*cosa + y*sina))/pow(1 + pow(tx*cosa - ty*sina,2),1.5); 

  // dDOCAW/dt0
  double t0shift = 4.0;  // ns
  double drift_shift = 0.0;
  if (drift_time < 0.0) {
    drift_shift = drift;
  } else {
    drift_shift =
        (du > 0.0 ? 1.0 : -1.0) *
        fdc_drift_distance(drift_time + t0shift);
  }
  derivatives[FDCTrackD::dW_dt0] = (drift_shift - drift) / t0shift;

// And the cathodes
//dDOCAW/ddeltax
derivatives[FDCTrackD::dDOCAC_dDeltaX] = 0.;

//dDOCAW/ddeltax
derivatives[FDCTrackD::dDOCAC_dDeltaZ] = ty*cosa + tx*sina;

//dDOCAW/ddeltaPhiX
derivatives[FDCTrackD::dDOCAC_dDeltaPhiX] = 0.;

//dDOCAW/ddeltaPhiX
derivatives[FDCTrackD::dDOCAC_dDeltaPhiY] = 0.;

//dDOCAW/ddeltaPhiX
derivatives[FDCTrackD::dDOCAC_dDeltaPhiZ] = -(x*cosa) + y*sina;

// dDOCAW/dx
derivatives[FDCTrackD::dDOCAC_dx] = sina;

// dDOCAW/dy
derivatives[FDCTrackD::dDOCAW_dy] = cosa;

// dDOCAW/dtx
derivatives[FDCTrackD::dDOCAW_dtx] = 0.;

// dDOCAW/dty
derivatives[FDCTrackD::dDOCAW_dty] = 0.;

thisPull.AddTrackDerivatives(derivatives);
    }
    pulls.push_back(thisPull);  
  }
  else if (best_trajectory[m].id>=1000){ // CDC Hit
    unsigned int id=best_trajectory[m].id-1000;
    A=cdc_updates[id].C*JT*C.Invert();
    Ss=cdc_updates[id].S+A*(Ss-S);
    
    dC=A*(Cs-C)*A.Transpose();
    Cs=cdc_updates[id].C+dC;
    if (VERBOSE > 10) {
jout << " In Smoothing Step Using ID " << id << "/" << cdc_updates.size() << " for ring " << cdchits[id]->wire->ring << endl;
jout << " A cdc_updates[id].C Ss Cs " << endl;
A.Print(); cdc_updates[id].C.Print(); Ss.Print(); Cs.Print();
    }
    if(!Cs.IsPosDef()) {
if (VERBOSE) jout << "Cs is not PosDef!" << endl;
return VALUE_OUT_OF_RANGE;
    }
    
    const DCDCWire *wire=cdchits[id]->wire;
    DVector3 origin=wire->origin;
    double z0=origin.z();
    double vz=wire->udir.z();
    DVector3 wdir=(1./vz)*wire->udir;
    DVector3 wirepos=origin+(best_trajectory[m].z-z0)*wdir;
    // Position and direction from state vector
    double x=Ss(state_x);
    double y=Ss(state_y);
    double tx=Ss(state_tx);
    double ty=Ss(state_ty);
    
    DVector3 pos0(x,y,best_trajectory[m].z);
    DVector3 tdir(tx,ty,1.);
    
    // Find the true doca to the wire
    DVector3 diff=pos0-origin;
    double dx0=diff.x(),dy0=diff.y();
    double wdir_dot_diff=diff.Dot(wdir);
    double tdir_dot_diff=diff.Dot(tdir);
    double tdir_dot_wdir=tdir.Dot(wdir);
    double tdir2=tdir.Mag2();
    double wdir2=wdir.Mag2();
    double D=tdir2*wdir2-tdir_dot_wdir*tdir_dot_wdir;
    double N=tdir_dot_wdir*wdir_dot_diff-wdir2*tdir_dot_diff;
    double N1=tdir2*wdir_dot_diff-tdir_dot_wdir*tdir_dot_diff;
    double scale=1./D;
    double s=scale*N;
    double t=scale*N1;
    diff+=s*tdir-t*wdir;
    double d=diff.Mag()+d_EPS; // prevent division by zero
    double ddrift = cdc_updates[id].ddrift;
    
    double resi = ddrift - d;
    
    
    // Track projection
    DMatrix1x4 H; DMatrix4x1 H_T;
    {
double one_over_d=1./d;
double diffx=diff.x(),diffy=diff.y(),diffz=diff.z();
double wx=wdir.x(),wy=wdir.y();

double dN1dtx=2.*tx*wdir_dot_diff-wx*tdir_dot_diff-tdir_dot_wdir*dx0;
double dDdtx=2.*tx*wdir2-2.*tdir_dot_wdir*wx;
double dtdtx=scale*(dN1dtx-t*dDdtx);

double dN1dty=2.*ty*wdir_dot_diff-wy*tdir_dot_diff-tdir_dot_wdir*dy0;
double dDdty=2.*ty*wdir2-2.*tdir_dot_wdir*wy;
double dtdty=scale*(dN1dty-t*dDdty);

double dNdtx=wx*wdir_dot_diff-wdir2*dx0;
double dsdtx=scale*(dNdtx-s*dDdtx);

double dNdty=wy*wdir_dot_diff-wdir2*dy0;
double dsdty=scale*(dNdty-s*dDdty);

H(state_tx)=H_T(state_tx)
 =one_over_d*(diffx*(s+tx*dsdtx-wx*dtdtx)+diffy*(ty*dsdtx-wy*dtdtx)
       +diffz*(dsdtx-dtdtx));
H(state_ty)=H_T(state_ty)
 =one_over_d*(diffx*(tx*dsdty-wx*dtdty)+diffy*(s+ty*dsdty-wy*dtdty)
       +diffz*(dsdty-dtdty));

double dsdx=scale*(tdir_dot_wdir*wx-wdir2*tx);
double dtdx=scale*(tdir2*wx-tdir_dot_wdir*tx);
double dsdy=scale*(tdir_dot_wdir*wy-wdir2*ty);
double dtdy=scale*(tdir2*wy-tdir_dot_wdir*ty);

H(state_x)=H_T(state_x)
 =one_over_d*(diffx*(1.+dsdx*tx-dtdx*wx)+diffy*(dsdx*ty-dtdx*wy)
       +diffz*(dsdx-dtdx));
H(state_y)=H_T(state_y)
 =one_over_d*(diffx*(dsdy*tx-dtdy*wx)+diffy*(1.+dsdy*ty-dtdy*wy)
       +diffz*(dsdy-dtdy));
    }
    double V=cdc_updates[id].V;

    if (VERBOSE > 10) jout << " d " << d << " H*S " << H*S << endl;
    if (cdchits[id]->wire->ring==RING_TO_SKIP){
V=V+H*Cs*H_T;
    }
    else{
V=V-H*Cs*H_T;
    }
    if (V<0) return VALUE_OUT_OF_RANGE;
    
    // Add the pull
    double myscale=1./sqrt(1.+tx*tx+ty*ty);
    double cosThetaRel=wire->udir.Dot(DVector3(myscale*tx,myscale*ty,myscale));
    DTrackFitter::pull_t thisPull(resi,sqrt(V),
		    best_trajectory[m].t*SPEED_OF_LIGHT29.9792458,
		    cdc_updates[id].tdrift,
		    d,cdchits[id], NULL__null,
		    diff.Phi(), //docaphi
		    best_trajectory[m].z,cosThetaRel,
		    cdc_updates[id].tdrift);

    // Derivatives for alignment
    double wtx=wire->udir.X(), wty=wire->udir.Y(), wtz=wire->udir.Z();
    double wx=wire->origin.X(), wy=wire->origin.Y(), wz=wire->origin.Z();
    
    double z=best_trajectory[m].z;
    double tx2=tx*tx, ty2=ty*ty;
    double wtx2=wtx*wtx, wty2=wty*wty, wtz2=wtz*wtz;
    double denom=(1 + ty2)*wtx2 + (1 + tx2)*wty2 - 2*ty*wty*wtz + (tx2 + ty2)*wtz2 - 2*tx*wtx*(ty*wty + wtz)+d_EPS;
    double denom2=denom*denom;
    double c1=-(wtx - tx*wtz)*(wy - y);
    double c2=wty*(wx - tx*wz - x + tx*z);
    double c3=ty*(-(wtz*wx) + wtx*wz + wtz*x - wtx*z);
    double dscale=0.5*(1/d);
    
    vector<double> derivatives(11);
    
    derivatives[CDCTrackD::dDOCAdOriginX]=dscale*(2*(wty - ty*wtz)*(c1 + c2 + c3))/denom;
    
    derivatives[CDCTrackD::dDOCAdOriginY]=dscale*(2*(-wtx + tx*wtz)*(c1 + c2 + c3))/denom;
    
    derivatives[CDCTrackD::dDOCAdOriginZ]=dscale*(2*(ty*wtx - tx*wty)*(c1 + c2 + c3))/denom;
    
    derivatives[CDCTrackD::dDOCAdDirX]=dscale*(2*(wty - ty*wtz)*(c1 + c2 + c3)*
				 (tx*(ty*wty + wtz)*(wx - x) + (wty - ty*wtz)*(-wy + y + ty*(wz - z)) +
				  wtx*(-((1 + ty2)*wx) + (1 + ty2)*x + tx*(ty*wy + wz - ty*y - z)) + tx2*(wty*(-wy + y) + wtz*(-wz + z))))/denom2;
    
    derivatives[CDCTrackD::dDOCAdDirY]=dscale*(-2*(wtx - tx*wtz)*(c1 + c2 + c3)*
				 (tx*(ty*wty + wtz)*(wx - x) + (wty - ty*wtz)*(-wy + y + ty*(wz - z)) +
				  wtx*(-((1 + ty2)*wx) + (1 + ty2)*x + tx*(ty*wy + wz - ty*y - z)) + tx2*(wty*(-wy + y) + wtz*(-wz + z))))/denom2;

    derivatives[CDCTrackD::dDOCAdDirZ]=dscale*(-2*(ty*wtx - tx*wty)*(c1 + c2 + c3)*
				 (-(tx*(ty*wty + wtz)*(wx - x)) + tx2*(wty*(wy - y) + wtz*(wz - z)) + (wty - ty*wtz)*(wy - y + ty*(-wz + z)) +
				  wtx*((1 + ty2)*wx - (1 + ty2)*x + tx*(-(ty*wy) - wz + ty*y + z))))/denom2;

    derivatives[CDCTrackD::dDOCAdS0]=-derivatives[CDCTrackD::dDOCAdOriginX];
    
    derivatives[CDCTrackD::dDOCAdS1]=-derivatives[CDCTrackD::dDOCAdOriginY];

    derivatives[CDCTrackD::dDOCAdS2]=dscale*(2*(wty - ty*wtz)*(-c1 - c2 - c3)*
			       (-(wtx*wtz*wx) - wty*wtz*wy + wtx2*wz + wty2*wz + wtx*wtz*x + wty*wtz*y - wtx2*z - wty2*z +
				tx*(wty2*(wx - x) + wtx*wty*(-wy + y) + wtz*(wtz*wx - wtx*wz - wtz*x + wtx*z)) +
				ty*(wtx*wty*(-wx + x) + wtx2*(wy - y) + wtz*(wtz*wy - wty*wz - wtz*y + wty*z))))/denom2;
    
    derivatives[CDCTrackD::dDOCAdS3]=dscale*(2*(wtx - tx*wtz)*(c1 + c2 + c3)*
			       (-(wtx*wtz*wx) - wty*wtz*wy + wtx2*wz + wty2*wz + wtx*wtz*x + wty*wtz*y - wtx2*z - wty2*z +
				tx*(wty2*(wx - x) + wtx*wty*(-wy + y) + wtz*(wtz*wx - wtx*wz - wtz*x + wtx*z)) +
				ty*(wtx*wty*(-wx + x) + wtx2*(wy - y) + wtz*(wtz*wy - wty*wz - wtz*y + wty*z))))/denom2;
    
    thisPull.AddTrackDerivatives(derivatives);
    
    pulls.push_back(thisPull);
  }
  else{
    A=best_trajectory[m].Ckk*JT*C.Invert();
    Ss=best_trajectory[m].Skk+A*(Ss-S);
    Cs=best_trajectory[m].Ckk+A*(Cs-C)*A.Transpose();
  }
  
  S=best_trajectory[m].Skk;
  C=best_trajectory[m].Ckk;
  JT=best_trajectory[m].J.Transpose();
}

return NOERROR;
2069}

2071shared_ptr<TMatrixFSym> 
2072DTrackFitterStraightTrack::Get7x7ErrorMatrix(TMatrixFSym C,DMatrix4x1 &S,double sign){
auto C7x7 = dResourcePool_TMatrixFSym->Get_SharedResource();
C7x7->ResizeTo(7, 7);
DMatrix J(7,5);

double p=5.; // fixed: cannot measure
double tx_=S(state_tx);
double ty_=S(state_ty);
double x_=S(state_x);
double y_=S(state_y);
double tanl=sign/sqrt(tx_*tx_+ty_*ty_);
double tanl2=tanl*tanl;
double lambda=atan(tanl);
double sinl=sin(lambda);
double sinl3=sinl*sinl*sinl;

J(state_X,state_x)=J(state_Y,state_y)=1.;
J(state_Pz,state_ty)=-p*ty_*sinl3;
J(state_Pz,state_tx)=-p*tx_*sinl3;
J(state_Px,state_ty)=J(state_Py,state_tx)=-p*tx_*ty_*sinl3;
J(state_Px,state_tx)=p*(1.+ty_*ty_)*sinl3;
J(state_Py,state_ty)=p*(1.+tx_*tx_)*sinl3;
J(state_Pz,4)=-p*p*sinl;
J(state_Px,4)=tx_*J(state_Pz,4);
J(state_Py,4)=ty_*J(state_Pz,4); 
J(state_Z,state_x)=-tx_*tanl2;
J(state_Z,state_y)=-ty_*tanl2;
double diff=tx_*tx_-ty_*ty_;
double frac=tanl2*tanl2;
J(state_Z,state_tx)=(x_*diff+2.*tx_*ty_*y_)*frac;
J(state_Z,state_ty)=(2.*tx_*ty_*x_-y_*diff)*frac;

// C'= JCJ^T
*C7x7=C.Similarity(J);

return C7x7;
2108}

2110// Routine to get extrapolations to other detectors
2111void DTrackFitterStraightTrack::GetExtrapolations(const DVector3 &pos0,
				  const DVector3 &dir){
double x0=pos0.x(),y0=pos0.y(),z0=pos0.z();
double s=0.,t=0.;
DVector3 pos(0,0,0);
DVector3 diff(0,0,0);
ClearExtrapolations();

double z=z0;
double dz=0.1;
double uz=dir.z();
double ux=dir.x()/uz,ux2=ux*ux;
double uy=dir.y()/uz,uy2=uy*uy;

// Extrapolate to start counter
double Rd=SC_BARREL_R;
double A=ux*x0 + uy*y0;
double B=uy2*(Rd - x0)*(Rd + x0) + 2*ux*uy*x0*y0 + ux2*(Rd - y0)*(Rd + y0);
double C=ux2 + uy2;
if (sc_pos.empty()==false && z<SC_END_NOSE_Z && B>0){
  dz=-(A-sqrt(B))/C;
  pos=pos0+(dz/uz)*dir;
  z=pos.z();
  bool done=(z<sc_pos[0][0].z()||z>SC_END_NOSE_Z);
  while(!done){
    double d_old=1000.,d=1000.;
    unsigned int index=0;
   
    for (unsigned int m=0;m<12;m++){
double dphi=pos.Phi()-SC_PHI_SECTOR1;
if (dphi<0) dphi+=2.*M_PI3.14159265358979323846;
index=int(floor(dphi/(2.*M_PI3.14159265358979323846/30.)));
if (index>29) index=0;
d=sc_norm[index][m].Dot(pos-sc_pos[index][m]);
if (d*d_old<0){ // break if we cross the current plane
 pos+=d*dir;
 s=(pos-pos0).Mag();
 t=s/29.98;
 extrapolations[SYS_START].push_back(DTrackFitter::Extrapolation_t(pos,dir,t,s));
 done=true;
 break;
}
d_old=d;
    }
    pos-=(0.1/uz)*dir;
    z=pos.z();
  }
}
// Extrapolate to BCAL
Rd=65.; // approximate BCAL inner radius
B=uy2*(Rd - x0)*(Rd + x0) + 2*ux*uy*x0*y0 + ux2*(Rd - y0)*(Rd + y0);
if (B>0){
  diff=-(A-sqrt(B))/(C*uz)*dir;
  s=diff.Mag();
  t=s/29.98;
  pos=pos0+diff;
  extrapolations[SYS_BCAL].push_back(DTrackFitter::Extrapolation_t(pos,dir,t,s)); 
  Rd=89.; // approximate BCAL outer radius
  B=uy2*(Rd - x0)*(Rd + x0) + 2*ux*uy*x0*y0 + ux2*(Rd - y0)*(Rd + y0);
  if (B>0){
    diff=-(A-sqrt(B))/(C*uz)*dir;
    s=diff.Mag();
    t=s/29.98;
    pos=pos0+diff;
    extrapolations[SYS_BCAL].push_back(DTrackFitter::Extrapolation_t(pos,dir,t,s));
  }
}

// Extrapolate to TRD
for (unsigned int i=0;i<dTRDz_vec.size();i++){
  diff=((dTRDz_vec[i]-z0)/uz)*dir;
  pos=pos0+diff;
  if (fabs(pos.x())<130. && fabs(pos.y())<130.){
    s=diff.Mag();
    t=s/29.98;
    extrapolations[SYS_TRD].push_back(DTrackFitter::Extrapolation_t(pos,dir,t,s));   	
  }
}

// Extrapolate to DIRC
diff=((dDIRCz-z0)/uz)*dir;
pos=pos0+diff;
if (fabs(pos.x())<130. && fabs(pos.y())<130.){
  s=diff.Mag();
  t=s/29.98;
  extrapolations[SYS_DIRC].push_back(DTrackFitter::Extrapolation_t(pos,dir,t,s));
}

// Extrapolate to TOF
diff=((dTOFz-z0)/uz)*dir;
pos=pos0+diff;
if (fabs(pos.x())<130. && fabs(pos.y())<130.){
  double s=diff.Mag();
  double t=s/29.98;
  s=diff.Mag();
  t=s/29.98;
  extrapolations[SYS_TOF].push_back(DTrackFitter::Extrapolation_t(pos,dir,t,s));	
}

// Extrapolate to FCAL
diff=((dFCALz-z0)/uz)*dir;
pos=pos0+diff;
if (fabs(pos.x())<130. && fabs(pos.y())<130.){
  s=diff.Mag();
  t=s/29.98;
  extrapolations[SYS_FCAL].push_back(DTrackFitter::Extrapolation_t(pos,dir,t,s));  

  // extrapolate to exit of FCAL
  diff=((dFCALz+45.-z0)/uz)*dir;
  pos=pos0+diff;
  s=diff.Mag();
  t=s/29.98;
  extrapolations[SYS_FCAL].push_back(DTrackFitter::Extrapolation_t(pos,dir,t,s));
} 
2225}

←

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

8#ifndef _JEventLoop_
9#define _JEventLoop_

11#include <sys/time.h>

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;

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>

37// The following is here just so we can use ROOT's THtml class to generate documentation.
38#include "cint.h"


41// Place everything in JANA namespace
42namespace jana{


45template<class T> class JFactory;
46class JApplication;
47class JEventProcessor;


50class JEventLoop{
public:

friend class JApplication;

enum data_source_t{
	DATA_NOT_AVAILABLE = 1,
	DATA_FROM_CACHE,
	DATA_FROM_SOURCE,
	DATA_FROM_FACTORY
};

typedef struct{
	string caller_name;
	string caller_tag;
	string callee_name;
	string callee_tag;
	double start_time;
	double end_time;
	data_source_t data_source;
}call_stack_t;

typedef struct{
	const char* factory_name;
	string tag;
	const char* filename;
	int line;
}error_call_stack_t;

                                  JEventLoop(JApplication *app); ///< Constructor
					   virtual ~JEventLoop(); ////< Destructor
		   virtual const char* className(void){return static_className();}
			static const char* static_className(void){return "JEventLoop";}

				 JApplication* GetJApplication(void) const {return app;} ///< Get pointer to the JApplication object
                                void RefreshProcessorListFromJApplication(void); ///< Re-copy the list of JEventProcessors from JApplication
                    virtual jerror_t AddFactory(JFactory_base* factory); ///< Add a factory
					  jerror_t RemoveFactory(JFactory_base* factory); ///< Remove a factory
                      JFactory_base* GetFactory(const string data_name, const char *tag="", bool allow_deftag=true); ///< Get a specific factory pointer
              vector<JFactory_base*> GetFactories(void) const {return factories;} ///< Get all factory pointers
                                void GetFactoryNames(vector<string> &factorynames); ///< Get names of all factories in name:tag format
                                void GetFactoryNames(map<string,string> &factorynames); ///< Get names of all factories in map with key=name, value=tag
                  map<string,string> GetDefaultTags(void) const {return default_tags;}
                            jerror_t ClearFactories(void); ///< Reset all factories in preparation for next event.
					  jerror_t PrintFactories(int sparsify=0); ///< Print a list of all factories.
                            jerror_t Print(const string data_name, const char *tag=""); ///< Print the data of the given type

				 JCalibration* GetJCalibration();
              template<class T> bool GetCalib(string namepath, map<string,T> &vals);
              template<class T> bool GetCalib(string namepath, vector<T> &vals);
              template<class T> bool GetCalib(string namepath, T &val);

                          JGeometry* GetJGeometry();
              template<class T> bool GetGeom(string namepath, map<string,T> &vals);
	    template<class T> bool GetGeom(string namepath, T &val);

                   JResourceManager* GetJResourceManager(void);
                              string GetResource(string namepath);
              template<class T> bool GetResource(string namepath, T vals, int event_number=0);

                                void Initialize(void); ///< Do initializations just before event processing starts
                            jerror_t Loop(void); ///< Loop over events
                            jerror_t OneEvent(uint64_t event_number); ///< Process a specific single event (if source supports it)
                            jerror_t OneEvent(void); ///< Process a single event
                         inline void Pause(void){pause = 1;} ///< Pause event processing
                         inline void Resume(void){pause = 0;} ///< Resume event processing
                         inline void Quit(void){quit = 1;} ///< Clean up and exit the event loop
                         inline bool GetQuit(void) const {return quit;}
                                void QuitProgram(void);

                               // Support for random access of events
                          bool HasRandomAccess(void);
                          void AddToEventQueue(uint64_t event_number){ next_events_to_process.push_back(event_number); }
                          void AddToEventQueue(list<uint64_t> &event_numbers) { next_events_to_process.insert(next_events_to_process.end(), event_numbers.begin(), event_numbers.end()); }
                list<uint64_t> GetEventQueue(void){ return next_events_to_process; }
                          void ClearEventQueue(void){ next_events_to_process.clear(); }

      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).
      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)
      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
          template<class T> jerror_t GetFromSource(vector<const T*> &t, JFactory_base *factory=NULL__null); ///< Get data object pointers from source.
                      inline JEvent& GetJEvent(void){return event;} ///< Get pointer to the current JEvent object.
                         inline void SetJEvent(JEvent *event){this->event = *event;} ///< Set the JEvent pointer.
                         inline void SetAutoFree(int auto_free){this->auto_free = auto_free;} ///< Set the Auto-Free flag on/off
                    inline pthread_t GetPThreadID(void) const {return pthread_id;} ///< Get the pthread of the thread to which this JEventLoop belongs
                              double GetInstantaneousRate(void) const {return rate_instantaneous;} ///< Get the current event processing rate
                              double GetIntegratedRate(void) const {return rate_integrated;} ///< Get the current event processing rate
                              double GetLastEventProcessingTime(void) const {return delta_time_single;}
                        unsigned int GetNevents(void) const {return Nevents;}

                         inline bool CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB);

                         inline bool GetCallStackRecordingStatus(void){ return record_call_stack; }
                         inline void DisableCallStackRecording(void){ record_call_stack = false; }
                         inline void EnableCallStackRecording(void){ record_call_stack = true; }
                     inline void CallStackStart(JEventLoop::call_stack_t &cs, const string &caller_name, const string &caller_tag, const string callee_name, const string callee_tag);
                     inline void CallStackEnd(JEventLoop::call_stack_t &cs);
         inline vector<call_stack_t> GetCallStack(void){return call_stack;} ///< Get the current factory call stack
                         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
                         inline void AddToErrorCallStack(error_call_stack_t &cs){error_call_stack.push_back(cs);} ///< Add layer to the factory call stack
   inline vector<error_call_stack_t> GetErrorCallStack(void){return error_call_stack;} ///< Get the current factory error call stack
                                void PrintErrorCallStack(void); ///< Print the current factory call stack

                      const JObject* FindByID(JObject::oid_t id); ///< Find a data object by its identifier.
          template<class T> const T* FindByID(JObject::oid_t id); ///< Find a data object by its type and identifier
                      JFactory_base* FindOwner(const JObject *t); ///< Find the factory that owns a data object by pointer
                      JFactory_base* FindOwner(JObject::oid_t id); ///< Find a factory that owns a data object by identifier

                                     // User defined references
              template<class T> void SetRef(T *t);    ///< Add a user reference to this JEventLoop (must be a pointer)
                template<class T> T* GetRef(void);   ///< Get a user-defined reference of a specific type
        template<class T> vector<T*> GetRefsT(void); ///< Get all user-defined refrences of a specicif type
   vector<pair<const char*, void*> > GetRefs(void){ return user_refs; }  ///< Get copy of full list of user-defined references
              template<class T> void RemoveRef(T *t); ///< Remove user reference from list

							   // Convenience methods wrapping JEvent methods of same name
                      uint64_t GetStatus(void){return event.GetStatus();}
                          bool GetStatusBit(uint32_t bit){return event.GetStatusBit(bit);}
                          bool SetStatusBit(uint32_t bit, bool val=true){return event.SetStatusBit(bit, val);}
                          bool ClearStatusBit(uint32_t bit){return event.ClearStatusBit(bit);}
                          void ClearStatus(void){event.ClearStatus();}
                          void SetStatusBitDescription(uint32_t bit, string description){event.SetStatusBitDescription(bit, description);}
                        string GetStatusBitDescription(uint32_t bit){return event.GetStatusBitDescription(bit);}
                          void GetStatusBitDescriptions(map<uint32_t, string> &status_bit_descriptions){return event.GetStatusBitDescriptions(status_bit_descriptions);}
                              string StatusWordToString(void);

private:
JEvent event;
vector<JFactory_base*> factories;
vector<JEventProcessor*> processors;
vector<error_call_stack_t> error_call_stack;
vector<call_stack_t> call_stack;
JApplication *app;
JThread *jthread;
bool initialized;
bool print_parameters_called;
int pause;
int quit;
int auto_free;
pthread_t pthread_id;
map<string, string> default_tags;
vector<pair<string,string> > auto_activated_factories;
bool record_call_stack;
string caller_name;
string caller_tag;
vector<uint64_t> event_boundaries;
int32_t event_boundaries_run; ///< Run number boundaries were retrieved from (possbily 0)
list<uint64_t> next_events_to_process;

uint64_t Nevents;			      ///< Total events processed (this thread)
uint64_t Nevents_rate;		   ///< Num. events accumulated for "instantaneous" rate
double delta_time_single;		///< Time spent processing last event
double delta_time_rate;			///< Integrated time accumulated "instantaneous" rate (partial number of events)
double delta_time;				///< Total time spent processing events (this thread)
double rate_instantaneous;		///< Latest instantaneous rate
double rate_integrated;			///< Rate integrated over all events
 
static data_source_t null_data_source;

vector<pair<const char*, void*> > user_refs;
210};


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

219#if !defined(__CINT__) && !defined(__CLING__)

221//-------------
222// GetSingle
223//-------------
224template<class T>
225JFactory<T>* JEventLoop::GetSingle(const T* &t, const char *tag, bool exception_if_not_one)
226{
/// This is a convenience method that can be used to get a pointer to the single
/// object of type T from the specified factory. It simply calls the Get(vector<...>) method
/// and copies the first pointer into "t" (or NULL if something other than 1 object is returned).
/// 
/// This is intended to address the common situation in which there is an interest
/// in the event if and only if there is exactly 1 object of type T. If the event
/// has no objects of that type or more than 1 object of that type (for the specified
/// factory) then an exception of type "unsigned long" is thrown with the value
/// being the number of objects of type T. You can supress the exception by setting
/// exception_if_not_one to false. In that case, you will have to check if t==NULL to
/// know if the call succeeded.
vector<const T*> v;
JFactory<T> *fac = Get(v, tag);

if(v.size()!=1){
t = NULL__null;
if(exception_if_not_one) throw v.size();
}

t = v[0];

return fac;
249}

251//-------------
252// Get
253//-------------
254template<class T> 
255JFactory<T>* JEventLoop::Get(vector<const T*> &t, const char *tag, bool allow_deftag)
256{
/// Retrieve or generate the array of objects of
/// type T for the curent event being processed
/// by this thread.
///
/// By default, preference is given to reading the
/// objects from the data source(e.g. file) before generating
/// them in the factory. A flag exists in the factory
/// however to change this so that the factory is
/// given preference.
///
/// Note that regardless of the setting of this flag,
/// the data are only either read in or generated once.
/// Ownership of the objects will always be with the
/// factory so subsequent calls will always return pointers to
/// the same data.
///
/// If the factory is called on to generate the data,
/// it is done by calling the factory's Get() method
/// which, in turn, calls the evnt() method. 
/// 
/// First, we just call the GetFromFactory() method.
/// It will make the initial decision as to whether
/// it should look in the source first or not. If
/// it returns NULL, then the factory couldn't be
/// found so we automatically try the file.
///
/// Note that if no factory exists to hold the objects
/// from the file, one can be created automatically
/// providing the <i>JANA:AUTOFACTORYCREATE</i>
/// configuration parameter is set.

// Check if a tag was specified for this data type to use for the
// default.
const char *mytag = tag1.1
'tag' is not equal to NULL
1
'tag' is not equal to NULL
1
'tag' is not equal to NULL
==NULL__null ? "":tag; // protection against NULL tags
2
←
'?' condition is false→
if(strlen(mytag)==0 && allow_deftag2.1
'allow_deftag' is true
1
'allow_deftag' is true
1
'allow_deftag' is true
){
3
←
Taking true branch→
map<string, string>::const_iterator iter = default_tags.find(T::static_className());
if(iter!=default_tags.end())tag = iter->second.c_str();
4
←
Assuming the condition is true→
5
←
Taking true branch→
6
←
Value assigned to 'tag'→
}


// If we are trying to keep track of the call stack then we
// need to add a new call_stack_t object to the the list
// and initialize it with the start time and caller/callee
// info.
call_stack_t cs;

// Optionally record starting info of call stack entry
if(record_call_stack) CallStackStart(cs, caller_name, caller_tag, T::static_className(), tag);
7
←
Assuming field 'record_call_stack' is false→
8
←
Taking false branch→

// Get the data (or at least try to)
JFactory<T>* factory=NULL__null;
try{
factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
9
←
Passing value via 2nd parameter 'tag'→
10
←
Calling 'JEventLoop::GetFromFactory'→
if(!factory){
	// No factory exists for this type and tag. It's possible
	// that the source may be able to provide the objects
	// but it will need a place to put them. We can create a
	// dumb JFactory just to hold the data in case the source
	// can provide the objects. Before we do though, make sure
	// the user condones this via the presence of the
	// "JANA:AUTOFACTORYCREATE" config parameter.
	string p;
	try{
		gPARMS->GetParameter("JANA:AUTOFACTORYCREATE", p);
	}catch(...){}
	if(p.size()==0){
		jout<<std::endl;
		_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;
		jout<<"No factory of type \""<<T::static_className()<<"\" with tag \""<<tag<<"\" exists."<<std::endl;
		jout<<"If you are reading objects from a file, I can auto-create a factory"<<std::endl;
		jout<<"of the appropriate type to hold the objects, but this feature is turned"<<std::endl;
		jout<<"off by default. To turn it on, set the \"JANA:AUTOFACTORYCREATE\""<<std::endl;
		jout<<"configuration parameter. This can usually be done by passing the"<<std::endl;
		jout<<"following argument to the program from the command line:"<<std::endl;
		jout<<std::endl;
		jout<<"   -PJANA:AUTOFACTORYCREATE=1"<<std::endl;
		jout<<std::endl;
		jout<<"Note that since the most commonly expected occurance of this situation."<<std::endl;
		jout<<"is an error, the program will now throw an exception so that the factory."<<std::endl;
		jout<<"call stack can be printed."<<std::endl;
		jout<<std::endl;
		throw exception();
	}else{
		AddFactory(new JFactory<T>(tag));
		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;
	
		// Now try once more. The GetFromFactory method will call
		// GetFromSource since it's empty.
		factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
	}
}
}catch(exception &e){
// Uh-oh, an exception was thrown. Add us to the call stack
// and re-throw the exception
error_call_stack_t ecs;
ecs.factory_name = T::static_className();
ecs.tag = tag;
ecs.filename = NULL__null;
error_call_stack.push_back(ecs);
throw e;
}

// If recording the call stack, update the end_time field and add to stack
if(record_call_stack) CallStackEnd(cs);

return factory;
363}

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{
// We need to find the factory providing data type T with
// tag given by "tag". 
vector<JFactory_base*>::iterator iter=factories.begin();
JFactory<T> *factory = NULL__null;
string className(T::static_className());

// Check if a tag was specified for this data type to use for the
// default.
const char *mytag = tag==NULL__null ? "":tag; // protection against NULL tags
11
←
Assuming 'tag' is equal to NULL→
12
←
Assuming pointer value is null→
13
←
'?' condition is true→
if(strlen(mytag)==0 && allow_deftag13.1
'allow_deftag' is true
1
'allow_deftag' is true
1
'allow_deftag' is true
){
14
←
Taking true branch→
map<string, string>::const_iterator iter = default_tags.find(className);
if(iter!=default_tags.end())tag = iter->second.c_str();
15
←
Assuming the condition is false→
16
←
Taking false branch→
}

for(; iter!=factories.end(); iter++){
17
←
Calling 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
20
←
Returning from 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
21
←
Loop condition is true.  Entering loop body→
// It turns out a long standing bug in g++ makes dynamic_cast return
// zero improperly when used on objects created on one side of
// a dynamically shared object (DSO) and the cast occurs on the 
// other side. I saw bug reports ranging from 2001 to 2004. I saw
// saw it first-hand on LinuxEL4 using g++ 3.4.5. This is too bad
// since it is much more elegant (and safe) to use dynamic_cast.
// To avoid this problem which can occur with plugins, we check
// the name of the data classes are the same. (sigh)
//factory = dynamic_cast<JFactory<T> *>(*iter);
if(className == (*iter)->GetDataClassName())factory = (JFactory<T>*)(*iter);
22
←
Taking true branch→
if(factory == NULL__null)continue;
23
←
Assuming 'factory' is not equal to NULL→
24
←
Taking false branch→
const char *factag = factory->Tag()==NULL__null ? "":factory->Tag();
25
←
Assuming the condition is true→
26
←
'?' condition is true→
if(!strcmp(factag, tag)){
27
←
Null pointer passed to 2nd parameter expecting 'nonnull'
	break;
}else{
	factory=NULL__null;
}
}

// If factory not found, just return now
if(!factory){
data_source = DATA_NOT_AVAILABLE;
return NULL__null;
}

// OK, we found the factory. If the evnt() routine has already
// been called, then just call the factory's Get() routine
// to return a copy of the existing data
if(factory->evnt_was_called()){
factory->CopyFrom(t);
data_source = DATA_FROM_CACHE;
return factory;
}

// Next option is to get the objects from the data source
if(factory->GetCheckSourceFirst()){
// If the object type/tag is found in the source, it
// will return NOERROR, even if there are zero instances
// of it. If it is not available in the source then it
// will return OBJECT_NOT_AVAILABLE.

jerror_t err = GetFromSource(t, factory);
if(err == NOERROR){
	// A return value of NOERROR means the source had the objects
	// even if there were zero of them.(If the source had no
	// information about the objects OBJECT_NOT_AVAILABLE would 
	// have been returned.)
	// The GetFromSource() call will eventually lead to a call to
	// the GetObjects() method of the concrete class derived
	// from JEventSource. That routine should copy the object
	// pointers into the factory using the factory's CopyTo()
	// method which also sets the evnt_called flag for the factory.
	// Note also that the "t" vector is then filled with a call
	// to the factory's CopyFrom() method in JEvent::GetObjects().
	// All we need to do now is just set the factory pointers in
	// the newly generated JObjects and return the factory pointer.

	factory->SetFactoryPointers();
	data_source = DATA_FROM_SOURCE;

	return factory;
}
}

// OK. It looks like we have to have the factory make this.
// Get pointers to data from the factory.
factory->Get(t);
factory->SetFactoryPointers();
data_source = DATA_FROM_FACTORY;

return factory;
457}

459//-------------
460// GetFromSource
461//-------------
462template<class T> 
463jerror_t JEventLoop::GetFromSource(vector<const T*> &t, JFactory_base *factory)
464{
/// This tries to get objects from the event source.
/// "factory" must be a valid pointer to a JFactory
/// object since that will take ownership of the objects
/// created by the source.
/// This should usually be called from JEventLoop::GetFromFactory
/// which is called from JEventLoop::Get. The latter will
/// create a dummy JFactory of the proper flavor and tag if
/// one does not already exist so if objects exist in the
/// file without a corresponding factory to create them, they
/// can still be used.
if(!factory)throw OBJECT_NOT_AVAILABLE;

return event.GetObjects(t, factory);
478}

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{
/// This is used to fill initial info into a call_stack_t stucture
/// for recording the call stack. It should be matched with a call
/// to CallStackEnd. It is normally called from the Get() method
/// above, but may also be used by external actors to manipulate
/// the call stack (presumably for good and not evil).

struct itimerval tmr;
getitimer(ITIMER_PROFITIMER_PROF, &tmr);

cs.caller_name    = this->caller_name;
cs.caller_tag     = this->caller_tag;
this->caller_name = cs.callee_name = callee_name;
this->caller_tag  = cs.callee_tag  = callee_tag;
cs.start_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
499}

501//-------------
502// CallStackEnd
503//-------------
504inline void JEventLoop::CallStackEnd(JEventLoop::call_stack_t &cs)
505{
/// Complete a call stack entry. This should be matched
/// with a previous call to CallStackStart which was
/// used to fill the cs structure.

struct itimerval tmr;
getitimer(ITIMER_PROFITIMER_PROF, &tmr);
cs.end_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
caller_name = cs.caller_name;
caller_tag  = cs.caller_tag;
call_stack.push_back(cs);
516}

518//-------------
519// CheckEventBoundary
520//-------------
521inline bool JEventLoop::CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB)
522{
/// Check whether the two event numbers span one or more boundaries
/// in the calibration/conditions database for the current run number.
/// Return true if they do and false if they don't. The first parameter
/// "event_numberA" is also checked if it lands on a boundary in which
/// case true is also returned. If event_numberB lands on a boundary,
/// but event_numberA does not, then false is returned.
///
/// This method is not expected to be called by a user. It is, however called,
/// everytime a JFactory's Get() method is called.

// Make sure our copy of the boundaries is up to date
if(event.GetRunNumber()!=event_boundaries_run){
event_boundaries.clear(); // in case we can't get the JCalibration pointer
JCalibration *jcalib = GetJCalibration();
if(jcalib)jcalib->GetEventBoundaries(event_boundaries);
event_boundaries_run = event.GetRunNumber();
}

// Loop over boundaries
for(unsigned int i=0; i<event_boundaries.size(); i++){
uint64_t eb = event_boundaries[i];
if((eb - event_numberA)*(eb - event_numberB) < 0.0 || eb==event_numberA){ // think about it ....
	// events span a boundary or is on a boundary. Return true
	return true;
}
}

return false;
551}

553//-------------
554// FindByID
555//-------------
556template<class T> 
557const T* JEventLoop::FindByID(JObject::oid_t id)
558{
/// This is a templated method that can be used in place
/// of the non-templated FindByID(oid_t) method if one knows
/// the class of the object with the specified id.
/// This method is faster than calling the non-templated
/// FindByID and dynamic_cast-ing the JObject since
/// this will only search the objects of factories that
/// produce the desired data type.
/// This method will cast the JObject pointer to one
/// of the specified type. To use this method,
/// a type is specified in the call as follows:
///
/// const DMyType *t = loop->FindByID<DMyType>(id);

// Loop over factories looking for ones that provide
// specified data type.
for(unsigned int i=0; i<factories.size(); i++){
if(factories[i]->GetDataClassName() != T::static_className())continue;

// This factory provides data of type T. Search it for
// the object with the specified id.
const JObject *my_obj = factories[i]->GetByID(id);
if(my_obj)return dynamic_cast<const T*>(my_obj);
}

return NULL__null;
584}

586//-------------
587// GetCalib (map)
588//-------------
589template<class T>
590bool JEventLoop::GetCalib(string namepath, map<string,T> &vals)
591{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

// Note that we could do this by making "vals" a generic type T thus, combining
// this with the vector version below. However, doing this explicitly will make
// it easier for the user to understand how to call us.

vals.clear();

JCalibration *calib = GetJCalibration();
if(!calib){
_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;
return true;
}

return calib->Get(namepath, vals, event.GetEventNumber());
608}

610//-------------
611// GetCalib (vector)
612//-------------
613template<class T> bool JEventLoop::GetCalib(string namepath, vector<T> &vals)
614{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

vals.clear();

JCalibration *calib = GetJCalibration();
if(!calib){
_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;
return true;
}

return calib->Get(namepath, vals, event.GetEventNumber());
627}

629//-------------
630// GetCalib (single)
631//-------------
632template<class T> bool JEventLoop::GetCalib(string namepath, T &val)
633{
/// This is a convenience method for getting a single entry. It
/// simply calls the vector version and returns the first entry.
/// It returns true if the vector version returns true AND there
/// is at least one entry in the vector. No check is made for there
/// there being more than one entry in the vector.

vector<T> vals;
bool ret = GetCalib(namepath, vals);
if(vals.empty()) return true;
val = vals[0];

return ret;
646}

648//-------------
649// GetGeom (map)
650//-------------
651template<class T>
652bool JEventLoop::GetGeom(string namepath, map<string,T> &vals)
653{
/// Get the JGeometry object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

// Note that we could do this by making "vals" a generic type T thus, combining
// this with the vector version below. However, doing this explicitly will make
// it easier for the user to understand how to call us.

vals.clear();

JGeometry *geom = GetJGeometry();
if(!geom){
_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;
return true;
}

return geom->Get(namepath, vals);
670}

672//-------------
673// GetGeom (atomic)
674//-------------
675template<class T> bool JEventLoop::GetGeom(string namepath, T &val)
676{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

JGeometry *geom = GetJGeometry();
if(!geom){
_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;
return true;
}

return geom->Get(namepath, val);
687}

689//-------------
690// SetRef
691//-------------
692template<class T>
693void JEventLoop::SetRef(T *t)
694{
pair<const char*, void*> p(typeid(T).name(), (void*)t);
user_refs.push_back(p);
697}

699//-------------
700// GetResource
701//-------------
702template<class T> bool JEventLoop::GetResource(string namepath, T vals, int event_number)
703{
JResourceManager *resource_manager = GetJResourceManager();
if(!resource_manager){
string mess = string("Unable to get the JResourceManager object (namepath=\"")+namepath+"\")";
throw JException(mess);
}

return resource_manager->Get(namepath, vals, event_number);
711}

713//-------------
714// GetRef
715//-------------
716template<class T>
717T* JEventLoop::GetRef(void)
718{
/// Get a user-defined reference (a pointer)
for(unsigned int i=0; i<user_refs.size(); i++){
if(user_refs[i].first == typeid(T).name()) return (T*)user_refs[i].second;
}

return NULL__null;
725}

727//-------------
728// GetRefsT
729//-------------
730template<class T>
731vector<T*> JEventLoop::GetRefsT(void)
732{
vector<T*> refs;
for(unsigned int i=0; i<user_refs.size(); i++){
if(user_refs[i].first == typeid(T).name()){
	refs.push_back((T*)user_refs[i].second);
}
}

return refs;
741}

743//-------------
744// RemoveRef
745//-------------
746template<class T>
747void JEventLoop::RemoveRef(T *t)
748{
vector<pair<const char*, void*> >::iterator iter;
for(iter=user_refs.begin(); iter!= user_refs.end(); iter++){
if(iter->second == (void*)t){
	user_refs.erase(iter);
	return;
}
}
_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}


760#endif //__CINT__  __CLING__

762} // Close JANA namespace



766#endif // _JEventLoop_


←

/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(); }
18
←
Assuming the condition is true→
19
←
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