TrkPocaXY Class Reference

#include <TrkPocaXY.h>

Inheritance diagram for TrkPocaXY:

List of all members.

Public Member Functions
	TrkPocaXY (const Trajectory &traj, double flt, const HepPoint3D &pt, double precision=1.0e-4)
	TrkPocaXY (const Trajectory &traj1, double flt1, const Trajectory &traj2, double flt2, double precision=1.0e-4)
	~TrkPocaXY ()
double	docaXY () const
double	fltl1 () const
double	fltl2 () const
const TrkErrCode &	status () const
double	flt1 () const
double	flt2 () const
double	precision ()
Protected Member Functions
void	minimize (const Trajectory &traj1, double f1, const Trajectory &traj2, double f2)
void	minimize (const Trajectory &traj1, double f1, const HepPoint3D &pt)
void	stepTowardPoca (const Trajectory &traj1, const Trajectory &traj2)
void	stepToPointPoca (const Trajectory &traj, const HepPoint3D &pt)
Protected Attributes
double	_precision
double	_flt1
double	_flt2
TrkErrCode	_status
Static Protected Attributes
static double	_maxDist = 1.e7
static int	_maxTry = 500
static double	_extrapToler = 2.
Private Member Functions
void	interLineCircle (const double &m, const double &q, const double &xc, const double &yc, const double &radius, double &xint1, double &yint1, double &xint2, double &yint2)
void	interTwoLines (const double &m1, const double &q1, const double &m2, const double &q2, double &xint, double &yint)
void	interTwoCircles (const double &xc1, const double &yc1, const double &r1, const double &xc2, const double &yc2, const double &r2, double &xint1, double &yint1, double &xint2, double &yint2)
Private Attributes
double	_docaxy

---

Detailed Description

Definition at line 27 of file TrkPocaXY.h.

---

Constructor & Destructor Documentation

TrkPocaXY::TrkPocaXY	(	const Trajectory &	traj,
		double	flt,
		const HepPoint3D &	pt,
		double	precision = 1.0e-4
	)

Definition at line 23 of file TrkPocaXY.cxx.

References _docaxy, TrkPocaBase::_flt1, TrkPocaBase::_flt2, TrkPocaBase::_status, and TrkPocaBase::status().

  : TrkPocaBase(flt,prec) , _docaxy(-9999.0) {
// construct a line trajectory through the given point
  Hep3Vector zaxis(0.0,0.0,1.0);
// length is set to 200cm (shouldn't matter)
//  TrkLineTraj zline(pt, zaxis, 200.0);//yzhang del
  TrkLineTraj zline(pt, zaxis, 2000.0);//yzhang change
// create a 2-traj poca with this
  double zlen(pt.z());
  TrkPoca zlinepoca(traj,flt,zline,zlen,prec);
// transfer over the data
  _status = zlinepoca.status();
  if(status().success()){
    _flt1 = zlinepoca.flt1();
    _flt2 = zlinepoca.flt2();
    _docaxy  = zlinepoca.doca(); // doca should be perpendicular to the zline so 2d by construction
  }
}

TrkPocaXY::TrkPocaXY	(	const Trajectory &	traj1,
		double	flt1,
		const Trajectory &	traj2,
		double	flt2,
		double	precision = 1.0e-4
	)

Definition at line 43 of file TrkPocaXY.cxx.

References _docaxy, TrkPocaBase::_flt1, TrkPocaBase::_flt2, TrkPocaBase::_status, Trajectory::curvature(), Trajectory::delDirect(), Trajectory::direction(), check_raw_filter::dist, TrkErrCode::failure(), TrkPocaBase::flt1(), interLineCircle(), interTwoCircles(), interTwoLines(), Trajectory::position(), q, TrkErrCode::setFailure(), TrkErrCode::setSuccess(), TrkPocaBase::status(), and TrkErrCode::success().

  : TrkPocaBase(flt1,flt2,prec) , _docaxy(-9999.0) {

  // this traj-traj version starts by approximating the trejectories as
  // either a line or a circle and treats separately the line-line, 
  // circle-circle, and line-circle cases.

  _flt1 = flt1;
  _flt2 = flt2;
  double delta1(10*prec);
  double delta2(10*prec);
  unsigned niter(0);
  static const unsigned maxiter(20);
  while(niter < maxiter && 
        ( delta1 > prec || delta2 > prec ) ){
    // get positions and directions
    HepPoint3D  pos1 = traj1.position(_flt1);
    HepPoint3D pos1b;
    if(niter == 0)  pos1b = pos1;
    HepPoint3D  pos2 = traj2.position(_flt2);
    Hep3Vector dir1 = traj1.direction(_flt1);
    Hep3Vector dir2 = traj2.direction(_flt2); 
    Hep3Vector dd1 = traj1.delDirect(_flt1);
    Hep3Vector dd2 = traj2.delDirect(_flt2);
    double curv1 = traj1.curvature(_flt1); 
    double curv2 = traj2.curvature(_flt2); 
    double m1, m2, q1, q2;
    double r1, r2, xc1,xc2,yc1,yc2,sr1,sr2;
    if(curv1 == 0) {
      // convert to m*x+q
      if(dir1[0] == 0){
        m1 = 1.0e12;
      }else{
        m1 = dir1[1]/dir1[0];
      }
      q1 = pos1.y()-pos1.x()*m1;
    }else{
    // get circle parameters
      r1 = (1- dir1[2]*dir1[2])/curv1;
      sr1=r1;
      if(dir1[0]*dd1[1] - dir1[1]*dd1[0] < 0) sr1 = -r1;
      double cosphi1 = dir1[0]/sqrt(dir1[0]*dir1[0]+dir1[1]*dir1[1]);
      double sinphi1 = cosphi1 * dir1[1]/dir1[0];
      xc1 = pos1.x() - sr1 * sinphi1;
      yc1 = pos1.y() + sr1 * cosphi1;
    }
    if(curv2 == 0) {
      if(dir2[0] == 0){
        m2 = 1.0e12;
      }else{ 
        m2 = dir2[1]/dir2[0];
      }
       q2 = pos2.y()-pos2.x()*m2;
    }else{
      r2 = (1-dir2[2]*dir2[2])/curv2;
      sr2 = r2;
      if(dir2[0]*dd2[1] - dir2[1]*dd2[0] < 0) sr2 = -r2;
      double cosphi2 = dir2[0]/sqrt(dir2[0]*dir2[0]+dir2[1]*dir2[1]);
      double sinphi2 = cosphi2 * dir2[1]/dir2[0];
      xc2 = pos2.x() - sr2 * sinphi2;
      yc2 = pos2.y() + sr2 * cosphi2;
    }
    double xint, yint, xint1, xint2, yint1, yint2, absdoca;
    _docaxy = 0;
    _status.setSuccess(3);

    // First the line-line case

    if(curv1==0 && curv2==0){
    // do the intersection in 2d
      interTwoLines(m1, q1, m2, q2, xint, yint);
    }

    //  next do the two circle case

    if(curv1 !=0 && curv2 !=0){
      //There are four cases
      double cdist = sqrt((xc1-xc2)*(xc1-xc2)+(yc1-yc2)*(yc1-yc2));
      // a - coincident centers
      if (fabs(cdist) < 1.e-12 ) {      
       //the algorithm fails because the points have all the same distance
        _status.setFailure(12, 
                  "TrkPocaXY:: the two circles have the same center...");
         return;
      }
     // b - Actual intersection
      if ( (fabs(r1-r2) < cdist) && (cdist < r1+r2 ) ) {      
        interTwoCircles(xc1,yc1,r1,xc2,yc2,r2,xint1,yint1,xint2,yint2);
        double dist1 = (xint1-pos1b.x())*(xint1-pos1b.x()) + 
                       (yint1-pos1b.y())*(yint1-pos1b.y());
        double dist2 = (xint2-pos1b.x())*(xint2-pos1b.x()) + 
                       (yint2-pos1b.y())*(yint2-pos1b.y());

        //choose closest to begining of track1
        if(dist1<dist2){
          xint = xint1;
          yint = yint1;
        } else {
          xint = xint2; 
          yint = yint2;
        }    
      }

    // c - nested circles and d - separated circles

      if ( (fabs(r1-r2) > cdist) ||  // nested circles
         (    cdist > (r1+r2) )) {   // separated circles
        // line going through the centers of the two circles

        double m = (yc1-yc2)/(xc1-xc2);  // y = m*x+q
        double q = yc1 - m*xc1;

        // intersection points between the line and the two circles

        double xint1, yint1, xint2, yint2, zOfDCrossT1;

        interLineCircle(m, q, xc1, yc1, r1, xint1, yint1, xint2, yint2);

        double xint3, yint3, xint4, yint4 ;
        interLineCircle(m, q, xc2, yc2, r2, xint3, yint3, xint4, yint4);
        if (fabs(r1-r2) > cdist) { // nested circles
          absdoca = fabs(r1-r2)-cdist;
#ifdef MDCPATREC_DEBUG
          cout << "ErrMsg(debugging) doing nested circles in TrkPocaXY " << endl;
#endif

          double dist1_3 = pow((xint1-xint3),2.) + pow((yint1-yint3),2.);
          double dist2_4 = pow((xint2-xint4),2.) + pow((yint2-yint4),2.);

          if(dist1_3<dist2_4) {
            xint = 0.5*(xint1+xint3);
            yint = 0.5*(yint1+yint3);
            zOfDCrossT1 = (xint3-xint1)*dir1[1]-(yint3-yint1)*dir1[0];
          } else {
            xint = 0.5*(xint2+xint4);
            yint = 0.5*(yint2+yint4);
            zOfDCrossT1 = (xint4-xint2)*dir1[1]-(yint4-yint2)*dir1[0];
          }  
          _docaxy = absdoca;
          if( zOfDCrossT1 > 0) _docaxy = -absdoca;   
        }
        if( cdist > (r1+r2) ) { //separated circles
          absdoca = cdist - (r1+r2);
#ifdef MDCPATREC_DEBUG
          cout << "ErrMsg(debugging) doing separated circles in TrkPocaXY"<< endl;
#endif

          double dist2_3 = pow((xint2-xint3),2.) + pow((yint2-yint3),2.);
          double dist1_4 = pow((xint1-xint4),2.) + pow((yint1-yint4),2.);
          if (dist2_3<dist1_4){
            xint = 0.5*(xint2+xint3);
            yint = 0.5*(yint2+yint3);
            zOfDCrossT1 = (xint3-xint1)*dir1[1]-(yint3-yint1)*dir1[0];
          } else {
            xint = 0.5*(xint1+xint4);
            yint = 0.5*(yint1+yint4);
            zOfDCrossT1 = (xint4-xint2)*dir1[1]-(yint4-yint2)*dir1[0];
          }
          _docaxy = absdoca;
          if( zOfDCrossT1 > 0) _docaxy = -absdoca;   
        }
      }
    }

    // Now do the line-circle case

    if((curv1 == 0 && curv2 !=0) || (curv1 != 0 && curv2 == 0)) {
  // distance between the line and the circle center
      HepVector dirT;
      double m,q,r,xc,yc, zOfDCrossT1;
      if(curv1==0) {m=m1; q=q1; r=r2; xc=xc2; yc=yc2; dirT=dir2;}
      else          {m=m2; q=q2; r=r1; xc=xc1; yc=yc1; dirT=dir1;}

      double dist = fabs((m*xc-yc+q)/sqrt(1+m*m));
      if(dist <= r) {

        // the intersection points

        interLineCircle(m,q,xc,yc,r, xint1, yint1, xint2, yint2);
        double dist1 = (xint1-pos1b.x())*(xint1-pos1b.x()) + 
                       (yint1-pos1b.y())*(yint1-pos1b.y());
        double dist2 = (xint2-pos1b.x())*(xint2-pos1b.x()) + 
                       (yint2-pos1b.y())*(yint2-pos1b.y());
        //choose closest to the beginning of track 1
        if(dist1<dist2){
          xint = xint1;
          yint = yint1;
        } else {
          xint = xint2;
          yint = yint2;
        }    
      } else { // no intersection points
#ifdef MDCPATREC_DEBUG
        cout << "ErrMsg(debugging) doing separated line-circle in TrkPocaXY"<<endl;
#endif

        // line going through the circle center and perpendicular to traj1

        double mperp = -1./m;
        double qperp = yc - mperp*xc;

        // intersection between this line and the two trajectories

        interLineCircle(mperp, qperp, xc, yc, r, xint1,yint1,xint2,yint2);

        double xint3,yint3;
        interTwoLines(m, q, mperp, qperp, xint3, yint3);

        double dist1_3 = pow((xint1-xint3),2.) + pow((yint1-yint3),2.);
        double dist2_3 = pow((xint2-xint3),2.) + pow((yint2-yint3),2.);
        if (dist1_3<dist2_3) {
          xint =  0.5*(xint3 + xint1);
          yint =  0.5*(yint3 + yint1);
          absdoca = sqrt((xint3-xint1)*(xint3-xint1)+
                         (yint3-yint1)*(yint3-yint1));
          zOfDCrossT1 = (xint3-xint1)*dirT[1]-(yint3-yint1)*dirT[0];
        }
        else {
          xint = 0.5*(xint3 + xint2);
          yint = 0.5*(yint3 + yint2);
          absdoca = sqrt((xint3-xint2)*(xint3-xint2)+
                         (yint3-yint2)*(yint3-yint2));
          zOfDCrossT1 = (xint3-xint2)*dirT[1]-(yint3-yint2)*dirT[0];
        }
         _docaxy = absdoca;
        if( zOfDCrossT1 > 0) _docaxy = -absdoca;   
      }
    }

    // get the flight lengths for the intersection point

    HepPoint3D intpt( xint, yint, 0.0);
    TrkPocaXY poca1(traj1,_flt1,intpt,prec);
    TrkPocaXY poca2(traj2,_flt2,intpt,prec);
    _status = poca2.status();
    if(poca1.status().success() && poca2.status().success()) {
       delta1 = fabs(_flt1 - poca1.flt1());
       delta2 = fabs(_flt2 - poca2.flt1());
       _flt1 = poca1.flt1();
       _flt2 = poca2.flt1(); 
    }else{
#ifdef MDCPATREC_WARNING
      cout << "ErrMsg(warning)" << " poca fialure " << endl;
#endif
      if(poca1.status().failure()) _status=poca1.status();
      break;
    }
    niter++;
  }
#ifdef MDCPATREC_DEBUG
  cout <<"ErrMsg(debugging)" <<"TrkPocaXY iterations = " << niter-1 
       << "       flight lengths " << _flt1 <<" " << _flt2 << endl 
       << "       deltas " << delta1 <<"  " << delta2 << endl;
#endif
  if(niter == maxiter){
#ifdef MDCPATREC_ROUTINE
    cout << "ErrMsg(routine)" << "TrkPocaXY:: did not converge" << endl;
#endif
    _status.setFailure(13, "TrkPocaXY:: did not converge"); 
  }
}

TrkPocaXY::~TrkPocaXY

(

)

[inline]

Definition at line 43 of file TrkPocaXY.h.

{}

---

Member Function Documentation

double TrkPocaXY::docaXY

(

)

const [inline]

Definition at line 73 of file TrkPocaXY.h.

References _docaxy.

{return _docaxy;}

double TrkPocaBase::flt1

(

)

const [inline, inherited]

Definition at line 65 of file TrkPocaBase.h.

References TrkPocaBase::_flt1.

Referenced by TrkDifPieceTraj::append(), TrkPoca::calcDist(), TrkDifPoca::calcDist(), TrkSimpTraj::changePoint(), fltl1(), TrkSimpleRep::getAllWeights(), TrkCompTrk::getAllWeights(), TrkPocaBase::minimize(), TrkDifPieceTraj::prepend(), TrkPocaBase::stepToPointPoca(), TrkPocaBase::stepTowardPoca(), TrkPoca::TrkPoca(), TrkPocaXY(), and TrkHitOnTrk::updatePoca().

{return _flt1;}

double TrkPocaBase::flt2

(

)

const [inline, inherited]

Definition at line 68 of file TrkPocaBase.h.

References TrkPocaBase::_flt2.

Referenced by TrkPoca::calcDist(), TrkDifPoca::calcDist(), MdcUtilitySvc::docaPatPar(), fltl2(), TrkPocaBase::minimize(), TrkPocaBase::stepTowardPoca(), and TrkHitOnTrk::updatePoca().

{return _flt2;}

double TrkPocaXY::fltl1

(

)

const [inline]

Definition at line 49 of file TrkPocaXY.h.

References TrkPocaBase::flt1().

{ return flt1(); }

double TrkPocaXY::fltl2

(

)

const [inline]

Definition at line 50 of file TrkPocaXY.h.

References TrkPocaBase::flt2().

{ return flt2(); }

void TrkPocaXY::interLineCircle	(	const double &	m,
		const double &	q,
		const double &	xc,
		const double &	yc,
		const double &	radius,
		double &	xint1,
		double &	yint1,
		double &	xint2,
		double &	yint2
	)			[private]

Definition at line 308 of file TrkPocaXY.cxx.

References TrkPocaBase::_status, alpha, and TrkErrCode::setFailure().

Referenced by TrkPocaXY().

{

 double alpha = 1 + m*m;

 double beta = -xc +m*(q-yc);

 double gamma = xc*xc + (q-yc)*(q-yc) -r*r;

 double Delta = beta*beta - alpha*gamma;

 if (Delta < 0) {

    _status.setFailure(14,
              "TrkPocaXY:: the line and the circle heve no intersections...");
    return;

  }
  else if (fabs(Delta) <1.e-12) {

  xint1 = -beta/alpha;
  xint2 = xint1;

  }
  else {

    double xPlus = -beta/alpha + sqrt(beta*beta - alpha*gamma)/alpha;
    double xMinus = -beta/alpha - sqrt(beta*beta - alpha*gamma)/alpha;

    if (xPlus > xMinus) {
      xint1 = xMinus;
      xint2 = xPlus;
    }
    else {
      xint1 = xPlus;
      xint2 = xMinus;
    }
  }
 yint1 = m*xint1 + q;
 yint2 = m*xint2 + q;

 return;
}

void TrkPocaXY::interTwoCircles	(	const double &	xc1,
		const double &	yc1,
		const double &	r1,
		const double &	xc2,
		const double &	yc2,
		const double &	r2,
		double &	xint1,
		double &	yint1,
		double &	xint2,
		double &	yint2
	)			[private]

Definition at line 358 of file TrkPocaXY.cxx.

References EvtCyclic3::A, alpha, EvtCyclic3::B, and EvtCyclic3::C.

Referenced by TrkPocaXY().

{

  double A = (xc1*xc1 + yc1*yc1 - r1*r1) - (xc2*xc2 + yc2*yc2 - r2*r2);

  double B = -xc1 + xc2;

  double C = -yc1 + yc2;

  double alpha = 1 + (B*B)/(C*C);

  double beta = -xc1 + B/C*(yc1+A/(2*C));

  double gamma = xc1*xc1 + (yc1+A/(2*C))*(yc1+A/(2*C)) - r1*r1;

  double Delta = beta*beta - alpha*gamma;

  if (Delta < 0) {

    _status.setFailure(14, "TrkPocaXY:: the two circles have no intersections..\
.");
    return;

  }
  else if (fabs(Delta) <1.e-12) {

    xint1 = -beta/alpha;
    xint2 = xint1;

  }
  else {
    double xPlus = -beta/alpha + sqrt(beta*beta - alpha*gamma)/alpha;
    double xMinus = -beta/alpha - sqrt(beta*beta - alpha*gamma)/alpha;

    if (xPlus > xMinus) {
      xint1 = xMinus;
      xint2 = xPlus;
    }
    else {
      xint1 = xPlus;
      xint2 = xMinus;
    }

  }

  yint1 = -(A+2*B*xint1)/(2*C);
  yint2 = -(A+2*B*xint2)/(2*C);


  return;
}

void TrkPocaXY::interTwoLines	(	const double &	m1,
		const double &	q1,
		const double &	m2,
		const double &	q2,
		double &	xint,
		double &	yint
	)			[private]

Definition at line 416 of file TrkPocaXY.cxx.

Referenced by TrkPocaXY().

{

  if (fabs(m1-m2) < 1.e-12) {  // parallel lines

    //the algorithm fails because the points have all the same distance

    _status.setFailure(13, "TrkPocaXY:: the two lines are parallel...");
    return;
  }
  else { // the lines have an intersection point

    xint = (q2-q1)/(m1-m2);
    yint = m1*xint + q1;
  }

return;
}

void TrkPocaBase::minimize	(	const Trajectory &	traj1,
		double	f1,
		const HepPoint3D &	pt
	)			[protected, inherited]

Definition at line 120 of file TrkPocaBase.cxx.

References TrkPocaBase::_flt1, TrkPocaBase::_maxTry, TrkPocaBase::_status, Trajectory::distTo1stError(), TrkPocaBase::flt1(), genRecEmupikp::i, Trajectory::position(), TrkPocaBase::precision(), TrkErrCode::setFailure(), TrkPocaBase::status(), TrkPocaBase::stepToPointPoca(), and TrkErrCode::succeed.

{
   _status=TrkErrCode::succeed;
  _flt1 = f1;
  int pathDir = 1;  // which way are we stepping (+/- 1)

  int nTinyStep = 0;  // number of consecutive tiny steps -- oscillation test
  int nOscills = 0;
  double fltLast = 0., fltBeforeLast = 0.; // another check for oscillation
  for (int i = 0; i < _maxTry; i++) {
    fltLast = flt1();
    stepToPointPoca(traj, pt);
    if (status().failure()) return;
    double step = flt1() - fltLast;
    pathDir = (step > 0.) ? 1 : -1;
    // Can we stop stepping?
    double distToErr = traj.distTo1stError(fltLast, precision(), pathDir);
    bool mustStep = (fabs(step) > distToErr && step != 0.);
    // Crude test for oscillation (around cusp point of piecewise traj, I hope)
    if (fabs(step) < 0.5 * precision()) {
      nTinyStep++;
    } else {
      nTinyStep = 0;
      if (i > 1) {
        if (fabs(step) >= fabs(fltBeforeLast-fltLast) &&
            fabs(fltBeforeLast-flt1()) <= fabs(step)) {
          nOscills++;
          double halfway = (fltBeforeLast + fltLast) / 2.;
          if ((traj.position(flt1()) - pt).mag() >
              (traj.position(halfway) - pt).mag()) _flt1 = halfway;
        }
      }
    }
    if (nTinyStep > 3) mustStep = false;
    if (nOscills > 2) {
      mustStep = false;
#ifdef MDCPATREC_WARNING
      std::cout<<" ErrMsg(warning) "<< "Alleged oscillation detected. "
              << step << "  " << fltLast-fltBeforeLast
              << "  " << i << "  " << std::endl;
#endif
    }
    if (!mustStep) return;
    fltBeforeLast = fltLast;
  }
  // Ran off the end of the loop
  _status.setFailure(2);
}

void TrkPocaBase::minimize	(	const Trajectory &	traj1,
		double	f1,
		const Trajectory &	traj2,
		double	f2
	)			[protected, inherited]

Definition at line 32 of file TrkPocaBase.cxx.

References TrkPocaBase::_flt1, TrkPocaBase::_flt2, TrkPocaBase::_maxTry, TrkPocaBase::_status, Trajectory::distTo1stError(), false, TrkPocaBase::flt1(), TrkPocaBase::flt2(), Trajectory::position(), TrkPocaBase::precision(), TrkErrCode::setFailure(), TrkErrCode::setSuccess(), TrkPocaBase::status(), TrkPocaBase::stepTowardPoca(), TrkErrCode::succeed, and true.

Referenced by TrkDifPoca::TrkDifPoca(), and TrkPoca::TrkPoca().

{
  // Last revision: Jan 2003, WDH
  const int maxnOscillStep = 5 ;
  const int maxnDivergingStep = 5 ;

  // initialize
  _status = TrkErrCode::succeed;
  _flt1 = f1;
  _flt2 = f2;
  
  static HepPoint3D newPos1, newPos2 ;
  double delta(0), prevdelta(0) ;
  int nOscillStep(0) ;
  int nDivergingStep(0) ;
  bool finished = false ;
  int istep(0) ;

  for (istep=0; istep < _maxTry && !finished; ++istep) {
    double prevflt1  = flt1();
    double prevflt2  = flt2() ;
    double prevprevdelta = prevdelta ;
    prevdelta = delta;
    
    stepTowardPoca(traj1, traj2);
    if( status().failure() ) {
      // failure in stepTowardPoca
      finished=true ;
    } else {
      newPos1 = traj1.position(flt1());
      newPos2 = traj2.position(flt2());
      delta = (newPos1 - newPos2).mag();
      double step1 = flt1() - prevflt1;
      double step2 = flt2() - prevflt2;
      int pathDir1 = (step1 > 0.) ? 1 : -1;
      int pathDir2 = (step2 > 0.) ? 1 : -1;
      
      // Can we stop stepping?
      double distToErr1 = traj1.distTo1stError(prevflt1, precision(), pathDir1);
      double distToErr2 = traj2.distTo1stError(prevflt2, precision(), pathDir2);
      
      // converged if very small steps, or if parallel
      finished = 
        (fabs(step1) < distToErr1 && fabs(step2) < distToErr2 ) ||
        (status().success() == 3) ;
 
      // we have to catch some problematic cases
      if( !finished && istep>2 && delta > prevdelta) {
        if( prevdelta > prevprevdelta) {
          // diverging
          if(++nDivergingStep>maxnDivergingStep) { 
            _status.setFailure(2) ; // code for `Failed to converge'
            finished = true ;
          }
        } else {
          nDivergingStep=0;
          // oscillating
          if(++nOscillStep>maxnOscillStep) {
            // bail out of oscillation. since the previous step was
            // better, use that one.
            _flt1 = prevflt1 ;
            _flt2 = prevflt2 ;
            _status.setSuccess(21, "Oscillating poca.") ;
            finished = true ;
          } else {
            // we might be oscillating, but we could also just have
            // stepped over the minimum. choose a solution `in
            // between'.
            _flt1 = prevflt1 + 0.5*step1 ;
            _flt2 = prevflt2 + 0.5*step2 ;          
            newPos1 = traj1.position(_flt1) ;
            newPos2 = traj2.position(_flt2) ;
            delta = (newPos1 - newPos2).mag() ;
          }
        }
      } 
    }
  }
  if(!finished) _status.setSuccess(2) ; // code for 'not converged' (yet)
}

double TrkPocaBase::precision

(

)

[inline, inherited]

Definition at line 59 of file TrkPocaBase.h.

References TrkPocaBase::_precision.

Referenced by TrkDifPoca::calcDist(), TrkPocaBase::minimize(), and TrkPocaBase::stepToPointPoca().

{return _precision;}

const TrkErrCode & TrkPocaBase::status

(

)

const [inline, inherited]

Definition at line 62 of file TrkPocaBase.h.

References TrkPocaBase::_status.

Referenced by TrkDifPieceTraj::append(), TrkPoca::calcDist(), TrkDifPoca::calcDist(), TrkSimpTraj::changePoint(), MdcHitOnTrack::dcaToWire(), TrkHitOnTrk::getFitStuff(), TrkPocaBase::minimize(), TrkDifPieceTraj::prepend(), TrkDifPoca::TrkDifPoca(), TrkPoca::TrkPoca(), TrkPocaXY(), and TrkHitOnTrk::updatePoca().

{return _status;}

void TrkPocaBase::stepToPointPoca	(	const Trajectory &	traj,
		const HepPoint3D &	pt
	)			[protected, inherited]

Definition at line 250 of file TrkPocaBase.cxx.

References TrkPocaBase::_extrapToler, TrkPocaBase::_flt1, TrkPocaBase::_maxDist, TrkPocaBase::_status, Trajectory::distTo2ndError(), TrkPocaBase::flt1(), Trajectory::getInfo(), TrkPocaBase::precision(), and TrkErrCode::setFailure().

Referenced by TrkPocaBase::minimize().

{
  // Unsightly uninitialized variables:
  static Hep3Vector dir, delDir;
  static HepPoint3D trajPos;

  traj.getInfo(flt1(), trajPos, dir, delDir);
  Hep3Vector delta = ((CLHEP::Hep3Vector)trajPos) - ((CLHEP::Hep3Vector)pt);
  double denom = 1. + delta.dot(delDir);
  if (fabs(denom)*_maxDist < 1. ) {
    _status.setFailure(11, "TrkPoca::ambiguous tight looper.");
    return;
  }
  double df = -delta.dot(dir) / fabs(denom);
  int pathDir = (df > 0.) ? 1 : -1;

  // Don't try going farther than worst parabolic approximation will allow:
  double distToErr = traj.distTo2ndError(flt1(), _extrapToler, pathDir);
  if (fabs(df)>distToErr) df = (df>0?distToErr:-distToErr);
  // Make the step slightly longer -- prevents quitting at kinks
  df += 0.001 * pathDir * precision();
  _flt1 += df;
}

void TrkPocaBase::stepTowardPoca	(	const Trajectory &	traj1,
		const Trajectory &	traj2
	)			[protected, inherited]

Definition at line 176 of file TrkPocaBase.cxx.

References TrkPocaBase::_extrapToler, TrkPocaBase::_flt1, TrkPocaBase::_flt2, TrkPocaBase::_maxDist, TrkPocaBase::_status, Trajectory::distTo2ndError(), TrkPocaBase::flt1(), TrkPocaBase::flt2(), Trajectory::getInfo(), and TrkErrCode::setSuccess().

Referenced by TrkPocaBase::minimize().

{ 
  // Last revision: Jan 2003, WDH
  
  // A bunch of unsightly uninitialized variables:
  static Hep3Vector dir1, dir2;
  static Hep3Vector delDir1, delDir2;
  static HepPoint3D pos1, pos2;

  traj1.getInfo(flt1(), pos1, dir1, delDir1);
  traj2.getInfo(flt2(), pos2, dir2, delDir2);
  Hep3Vector delta = ((CLHEP::Hep3Vector)pos1) - ((CLHEP::Hep3Vector)pos2);
  double ua = -delta.dot(dir1);
  double ub =  delta.dot(dir2);
  double caa = dir1.mag2() + delta.dot(delDir1);
  double cbb = dir2.mag2() - delta.dot(delDir2);
  double cab = -dir1.dot(dir2);
  double det = caa * cbb - cab * cab;

  if(det<0) {
    // get rid of second order terms
    caa = dir1.mag2() ;
    cbb = dir2.mag2() ;
    det = caa * cbb - cab * cab;
  }
          
  if ( det < 1.e-8) {
    // If they are parallel (in quadratic approximation) give up
    _status.setSuccess(3);
    return;
  }

  double df1 = (ua * cbb - ub * cab)/det;
  int pathDir1 = (df1 > 0) ? 1 : -1;
  double df2 = (ub * caa - ua * cab)/det;
  int pathDir2 = (df2 > 0) ? 1 : -1;

  // Don't try going farther than worst parabolic approximation will
  // allow: Since ` _extrapToler' is large, this cut effectively only
  // takes care that we don't make large jumps past the kink in a
  // piecewise trajectory.

  double distToErr1 = traj1.distTo2ndError(flt1(), _extrapToler, pathDir1);
  double distToErr2 = traj2.distTo2ndError(flt2(), _extrapToler, pathDir2);

  // Factor to push just over border of piecewise traj (essential!)
  const double smudge = 1.01 ; 
  if( fabs(df1) > smudge*distToErr1 ) {
    // choose solution for which df1 steps just over border
    df1 = smudge*distToErr1 * pathDir1 ;
    // now recalculate df2, given df1:
    df2 = (ub - df1*cab)/cbb ;
  }

  if( fabs(df2) > smudge*distToErr2 ) {
    // choose solution for which df2 steps just over border
    df2 = smudge*distToErr2 * pathDir2 ;
    // now recalculate df1, given df2:
    df1 = (ua - df2*cab)/cbb ;
    // if still not okay,
    if( fabs(df1) > smudge*distToErr1 ) {
      df1 = smudge*distToErr1 * pathDir1 ;
    }
  }

  _flt1 += df1;
  _flt2 += df2;

  // another check for parallel trajectories
  if (fabs(flt1()) > _maxDist && fabs(flt2()) > _maxDist) 
    _status.setSuccess(3) ;
}

---

Member Data Documentation

double TrkPocaXY::_docaxy [private]

Definition at line 54 of file TrkPocaXY.h.

Referenced by docaXY(), and TrkPocaXY().

double TrkPocaBase::_extrapToler = 2. [static, protected, inherited]

Definition at line 54 of file TrkPocaBase.h.

Referenced by TrkPocaBase::stepToPointPoca(), and TrkPocaBase::stepTowardPoca().

double TrkPocaBase::_flt1 [protected, inherited]

Definition at line 40 of file TrkPocaBase.h.

Referenced by TrkPocaBase::flt1(), TrkPocaBase::minimize(), TrkPocaBase::stepToPointPoca(), TrkPocaBase::stepTowardPoca(), and TrkPocaXY().

double TrkPocaBase::_flt2 [protected, inherited]

Definition at line 41 of file TrkPocaBase.h.

Referenced by TrkPocaBase::flt2(), TrkPocaBase::minimize(), TrkPocaBase::stepTowardPoca(), and TrkPocaXY().

double TrkPocaBase::_maxDist = 1.e7 [static, protected, inherited]

Definition at line 52 of file TrkPocaBase.h.

Referenced by TrkPocaBase::stepToPointPoca(), and TrkPocaBase::stepTowardPoca().

int TrkPocaBase::_maxTry = 500 [static, protected, inherited]

Definition at line 53 of file TrkPocaBase.h.

Referenced by TrkPocaBase::minimize().

double TrkPocaBase::_precision [protected, inherited]

Definition at line 39 of file TrkPocaBase.h.

Referenced by TrkPocaBase::precision().

TrkErrCode TrkPocaBase::_status [protected, inherited]

Definition at line 42 of file TrkPocaBase.h.

Referenced by TrkDifPoca::calcDist(), interLineCircle(), TrkPocaBase::minimize(), TrkPocaBase::status(), TrkPocaBase::stepToPointPoca(), TrkPocaBase::stepTowardPoca(), and TrkPocaXY().

---