TrkDifPoca Class Reference

#include <TrkDifPoca.h>

Inheritance diagram for TrkDifPoca:

List of all members.

Public Member Functions
	TrkDifPoca (const TrkDifTraj &traj1, double flt1, const Trajectory &traj2, double flt2, double precision=1.e-5)
	TrkDifPoca (const TrkDifTraj &traj, double flt, const HepPoint3D &pt, double precision=1.e-5)
	~TrkDifPoca ()
const DifNumber &	difDoca () const
double	doca () const
const HepVector	derivs () const
void	fetchDerivs (HepVector &) 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	calcDist (const TrkDifTraj &traj1, const Trajectory &traj2)
void	calcDist (const TrkDifTraj &, const HepPoint3D &)
Private Attributes
DifNumber	_doca

---

Detailed Description

Definition at line 31 of file TrkDifPoca.h.

---

Constructor & Destructor Documentation

TrkDifPoca::TrkDifPoca	(	const TrkDifTraj &	traj1,
		double	flt1,
		const Trajectory &	traj2,
		double	flt2,
		double	precision = 1.e-5
	)

Definition at line 25 of file TrkDifPoca.cxx.

References calcDist(), TrkPocaBase::minimize(), and TrkPocaBase::status().

        : TrkPocaBase(f1,f2,prec), _doca(-9999.,0) 
{
        minimize(traj1,f1,traj2,f2);
        if (status().failure()) return;
        calcDist(traj1,traj2);

}

TrkDifPoca::TrkDifPoca	(	const TrkDifTraj &	traj,
		double	flt,
		const HepPoint3D &	pt,
		double	precision = 1.e-5
	)

Definition at line 36 of file TrkDifPoca.cxx.

References calcDist(), TrkPocaBase::minimize(), and TrkPocaBase::status().

        : TrkPocaBase(f1,prec), _doca(-9999.,0) 
{ 
        minimize(traj,f1,pt);
        if (status().failure()) return;
        calcDist(traj,pt);
}

TrkDifPoca::~TrkDifPoca

(

)

[inline]

Definition at line 39 of file TrkDifPoca.h.

{};

---

Member Function Documentation

void TrkDifPoca::calcDist	(	const TrkDifTraj &	,
		const HepPoint3D &
	)			[private]

Definition at line 85 of file TrkDifPoca.cxx.

References _doca, TrkPocaBase::_status, check_raw_filter::dist, TrkPocaBase::flt1(), TrkDifTraj::getDFInfo2(), DifVector::length(), DifNumber::number(), TrkPocaBase::precision(), TrkErrCode::setSuccess(), DifVector::x, and DifVector::z.

{
  // Does a final calculation of the distance -- and handles singularity 
  // in derivs when d = 0.

  // A bunch of unsightly uninitialized variables:
  static DifVector dir;
  static DifPoint posTraj;

  DifPoint posPoint(pt);
  traj.getDFInfo2(flt1(), posTraj, dir);  // call faster one, if exists
  DifVector delta = posTraj - posPoint;
  // delta -= dir*delta;  // note: dir*delta is zero, but the derivatives are NOT

  DifNumber dist = delta.length();
  if (dist.number() > 0.01 * precision()) { // d != 0
    _doca = dist;
  } else {
    // d = 0.  Fudge like mad: pick a direction (not parallel to traj) and 
    //   move the point slightly.  This should not happen very often.
    Hep3Vector fudgeVec(0., 0.1 * precision(), 0.0);
    if (dir.x.number() == 0.0 && dir.z.number() == 0.0) {
      fudgeVec = Hep3Vector(0.1 * precision(), 0., 0.);
    }
    posPoint += fudgeVec;
    delta = posTraj - posPoint;
    _doca = dist;
    _status.setSuccess(20, "TrkDifPoca: distance zero, calculation fudged.");
  }
}

void TrkDifPoca::calcDist	(	const TrkDifTraj &	traj1,
		const Trajectory &	traj2
	)			[private]

Definition at line 47 of file TrkDifPoca.cxx.

References _doca, TrkPocaBase::_status, cross(), DifVector::dot(), DifNumber::flipsign(), TrkPocaBase::flt1(), TrkPocaBase::flt2(), TrkDifTraj::getDFInfo2(), Trajectory::getInfo(), DifVector::normalize(), DifNumber::number(), TrkPocaBase::precision(), TrkErrCode::setFailure(), and TrkPocaBase::status().

Referenced by TrkDifPoca().

{
  // Does a final calculation of the distance -- better behaved derivs than 
  //  stepTowardPoca for zero dist.  In case of (nearly) parallel, returns 
  //  distance calculated at whatever point we happen to be at.  
  //  Derivatives are taken (of dist) w/r/t traj2 
  //  parameters, using DifNumbers for the calculation.

  // A bunch of unsightly uninitialized variables:
  static DifVector dir1;
  static DifPoint pos1;
  static Hep3Vector dirTem2;
  static HepPoint3D posTem2;

  // Request DifNumber info from traj1, ordinary info from traj2, and then 
  //   convert ordinary info into DifNumber
  traj2.getInfo(flt2(), posTem2, dirTem2);
  traj1.getDFInfo2(flt1(), pos1, dir1);
  DifVector dir2(dirTem2);
  DifPoint pos2(posTem2);

  DifVector delta = pos2 - pos1;
  if (status().success() != 3) { // Not parallel:
    DifVector between = cross( dir1, dir2 );  // cross-product
    between.normalize();
    _doca = delta * between;
  } else {  // Parallel: Arbitrary sign convention
    _doca = (delta - delta.dot(dir1) * dir1).length();
    if  (dir1.dot(dir2) > 0.) _doca.flipsign();
    if (fabs(_doca.number()) < 0.0001 * precision()) {
      // Parallel and on top of each other (derivatives singular)
      _doca = 0;
      _status.setFailure(3);
    }
  }
}

const HepVector TrkDifPoca::derivs

(

)

const [inline]

Definition at line 66 of file TrkDifPoca.h.

References _doca, and DifNumber::derivatives().

Referenced by TrkHitOnTrk::getFitStuff().

{return _doca.derivatives();}

const DifNumber & TrkDifPoca::difDoca

(

)

const [inline]

Definition at line 65 of file TrkDifPoca.h.

References _doca.

{return _doca;}

double TrkDifPoca::doca

(

)

const [inline]

Definition at line 64 of file TrkDifPoca.h.

References _doca, and DifNumber::number().

{return _doca.number();}

void TrkDifPoca::fetchDerivs

(

HepVector &

)

const [inline]

Definition at line 67 of file TrkDifPoca.h.

References _doca, and DifNumber::fetchDerivatives().

Referenced by TrkHitOnTrk::getFitStuff().

{_doca.fetchDerivatives(dv);}

double TrkPocaBase::flt1

(

)

const [inline, inherited]

Definition at line 65 of file TrkPocaBase.h.

References TrkPocaBase::_flt1.

Referenced by TrkDifPieceTraj::append(), TrkPoca::calcDist(), calcDist(), TrkSimpTraj::changePoint(), TrkPocaXY::fltl1(), TrkSimpleRep::getAllWeights(), TrkCompTrk::getAllWeights(), TrkPocaBase::minimize(), TrkDifPieceTraj::prepend(), TrkPocaBase::stepToPointPoca(), TrkPocaBase::stepTowardPoca(), TrkPoca::TrkPoca(), TrkPocaXY::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(), calcDist(), MdcUtilitySvc::docaPatPar(), TrkPocaXY::fltl2(), TrkPocaBase::minimize(), TrkPocaBase::stepTowardPoca(), and TrkHitOnTrk::updatePoca().

{return _flt2;}

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(), 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 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(), calcDist(), TrkSimpTraj::changePoint(), MdcHitOnTrack::dcaToWire(), TrkHitOnTrk::getFitStuff(), TrkPocaBase::minimize(), TrkDifPieceTraj::prepend(), TrkDifPoca(), TrkPoca::TrkPoca(), TrkPocaXY::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

DifNumber TrkDifPoca::_doca [private]

Definition at line 53 of file TrkDifPoca.h.

Referenced by calcDist(), derivs(), difDoca(), doca(), and fetchDerivs().

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::TrkPocaXY().

double TrkPocaBase::_flt2 [protected, inherited]

Definition at line 41 of file TrkPocaBase.h.

Referenced by TrkPocaBase::flt2(), TrkPocaBase::minimize(), TrkPocaBase::stepTowardPoca(), and TrkPocaXY::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 calcDist(), TrkPocaXY::interLineCircle(), TrkPocaBase::minimize(), TrkPocaBase::status(), TrkPocaBase::stepToPointPoca(), TrkPocaBase::stepTowardPoca(), and TrkPocaXY::TrkPocaXY().

---