Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

TrkDifPoca Class Reference

#include <TrkDifPoca.h>

Inheritance diagram for TrkDifPoca:

TrkPocaBase TrkPocaBase List of all members.

Public Member Functions

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

Protected Member Functions

void minimize (const Trajectory &traj1, double f1, const HepPoint3D &pt)
void minimize (const Trajectory &traj1, double f1, const Trajectory &traj2, double f2)
void minimize (const Trajectory &traj1, double f1, const HepPoint3D &pt)
void minimize (const Trajectory &traj1, double f1, const Trajectory &traj2, double f2)
void stepToPointPoca (const Trajectory &traj, const HepPoint3D &pt)
void stepToPointPoca (const Trajectory &traj, const HepPoint3D &pt)
void stepTowardPoca (const Trajectory &traj1, const Trajectory &traj2)
void stepTowardPoca (const Trajectory &traj1, const Trajectory &traj2)

Protected Attributes

double _flt1
double _flt2
double _precision
TrkErrCode _status

Static Protected Attributes

double _extrapToler = 2.
double _maxDist = 1.e7
int _maxTry = 500

Private Member Functions

void calcDist (const TrkDifTraj &, const HepPoint3D &)
void calcDist (const TrkDifTraj &traj1, const Trajectory &traj2)
void calcDist (const TrkDifTraj &, const HepPoint3D &)
void calcDist (const TrkDifTraj &traj1, const Trajectory &traj2)

Private Attributes

DifNumber _doca

Constructor & Destructor Documentation

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

00027         : TrkPocaBase(f1,f2,prec), _doca(-9999.,0) 
00028 {
00029         minimize(traj1,f1,traj2,f2);
00030         if (status().failure()) return;
00031         calcDist(traj1,traj2);
00032 
00033 }

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

00038         : TrkPocaBase(f1,prec), _doca(-9999.,0) 
00039 { 
00040         minimize(traj,f1,pt);
00041         if (status().failure()) return;
00042         calcDist(traj,pt);
00043 }

TrkDifPoca::~TrkDifPoca (   )  [inline]
  

00039 {};

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

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

TrkDifPoca::~TrkDifPoca (   )  [inline]
  

00039 {};

Member Function Documentation

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

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

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

00086 {
00087   // Does a final calculation of the distance -- and handles singularity 
00088   // in derivs when d = 0.
00089 
00090   // A bunch of unsightly uninitialized variables:
00091   static DifVector dir;
00092   static DifPoint posTraj;
00093 
00094   DifPoint posPoint(pt);
00095   traj.getDFInfo2(flt1(), posTraj, dir);  // call faster one, if exists
00096   DifVector delta = posTraj - posPoint;
00097   // delta -= dir*delta;  // note: dir*delta is zero, but the derivatives are NOT
00098 
00099   DifNumber dist = delta.length();
00100   if (dist.number() > 0.01 * precision()) { // d != 0
00101     _doca = dist;
00102   } else {
00103     // d = 0.  Fudge like mad: pick a direction (not parallel to traj) and 
00104     //   move the point slightly.  This should not happen very often.
00105     Hep3Vector fudgeVec(0., 0.1 * precision(), 0.0);
00106     if (dir.x.number() == 0.0 && dir.z.number() == 0.0) {
00107       fudgeVec = Hep3Vector(0.1 * precision(), 0., 0.);
00108     }
00109     posPoint += fudgeVec;
00110     delta = posTraj - posPoint;
00111     _doca = dist;
00112     _status.setSuccess(20, "TrkDifPoca: distance zero, calculation fudged.");
00113   }
00114 }

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

00048 {
00049   // Does a final calculation of the distance -- better behaved derivs than 
00050   //  stepTowardPoca for zero dist.  In case of (nearly) parallel, returns 
00051   //  distance calculated at whatever point we happen to be at.  
00052   //  Derivatives are taken (of dist) w/r/t traj2 
00053   //  parameters, using DifNumbers for the calculation.
00054 
00055   // A bunch of unsightly uninitialized variables:
00056   static DifVector dir1;
00057   static DifPoint pos1;
00058   static Hep3Vector dirTem2;
00059   static HepPoint3D posTem2;
00060 
00061   // Request DifNumber info from traj1, ordinary info from traj2, and then 
00062   //   convert ordinary info into DifNumber
00063   traj2.getInfo(flt2(), posTem2, dirTem2);
00064   traj1.getDFInfo2(flt1(), pos1, dir1);
00065   DifVector dir2(dirTem2);
00066   DifPoint pos2(posTem2);
00067 
00068   DifVector delta = pos2 - pos1;
00069   if (status().success() != 3) { // Not parallel:
00070     DifVector between = cross( dir1, dir2 );  // cross-product
00071     between.normalize();
00072     _doca = delta * between;
00073   } else {  // Parallel: Arbitrary sign convention
00074     _doca = (delta - delta.dot(dir1) * dir1).length();
00075     if  (dir1.dot(dir2) > 0.) _doca.flipsign();
00076     if (fabs(_doca.number()) < 0.0001 * precision()) {
00077       // Parallel and on top of each other (derivatives singular)
00078       _doca = 0;
00079       _status.setFailure(3);
00080     }
00081   }
00082 }

const HepVector TrkDifPoca::derivs (   )  const [inline]
 

const HepVector TrkDifPoca::derivs (   )  const [inline]
  

00066 {return _doca.derivatives();}

const DifNumber& TrkDifPoca::difDoca (   )  const [inline]
 

const DifNumber & TrkDifPoca::difDoca (   )  const [inline]
  

00065 {return _doca;}

double TrkDifPoca::doca (   )  const [inline]
 

double TrkDifPoca::doca (   )  const [inline]
  

00064 {return _doca.number();}

void TrkDifPoca::fetchDerivs (  HepVector &   )  const [inline]
 

void TrkDifPoca::fetchDerivs (  HepVector &   )  const [inline]
  

00067 {_doca.fetchDerivatives(dv);}

double TrkPocaBase::flt1 (   )  const [inline, inherited]
 

double TrkPocaBase::flt1 (   )  const [inline, inherited]
  

00065 {return _flt1;}

double TrkPocaBase::flt2 (   )  const [inline, inherited]
 

double TrkPocaBase::flt2 (   )  const [inline, inherited]
  

00068 {return _flt2;}

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

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

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

00122 {
00123    _status=TrkErrCode::succeed;
00124   _flt1 = f1;
00125   int pathDir = 1;  // which way are we stepping (+/- 1)
00126 
00127   int nTinyStep = 0;  // number of consecutive tiny steps -- oscillation test
00128   int nOscills = 0;
00129   double fltLast = 0., fltBeforeLast = 0.; // another check for oscillation
00130   for (int i = 0; i < _maxTry; i++) {
00131     fltLast = flt1();
00132     stepToPointPoca(traj, pt);
00133     if (status().failure()) return;
00134     double step = flt1() - fltLast;
00135     pathDir = (step > 0.) ? 1 : -1;
00136     // Can we stop stepping?
00137     double distToErr = traj.distTo1stError(fltLast, precision(), pathDir);
00138     bool mustStep = (fabs(step) > distToErr && step != 0.);
00139     // Crude test for oscillation (around cusp point of piecewise traj, I hope)
00140     if (fabs(step) < 0.5 * precision()) {
00141       nTinyStep++;
00142     } else {
00143       nTinyStep = 0;
00144       if (i > 1) {
00145         if (fabs(step) >= fabs(fltBeforeLast-fltLast) &&
00146             fabs(fltBeforeLast-flt1()) <= fabs(step)) {
00147           nOscills++;
00148           double halfway = (fltBeforeLast + fltLast) / 2.;
00149           if ((traj.position(flt1()) - pt).mag() >
00150               (traj.position(halfway) - pt).mag()) _flt1 = halfway;
00151         }
00152       }
00153     }
00154     if (nTinyStep > 3) mustStep = false;
00155     if (nOscills > 2) {
00156       mustStep = false;
00157 #ifdef MDCPATREC_WARNING
00158       std::cout<<" ErrMsg(warning) "<< "Alleged oscillation detected. "
00159               << step << "  " << fltLast-fltBeforeLast
00160               << "  " << i << "  " << std::endl;
00161 #endif
00162     }
00163     if (!mustStep) return;
00164     fltBeforeLast = fltLast;
00165   }
00166   // Ran off the end of the loop
00167   _status.setFailure(2);
00168 }

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

00034 {
00035   // Last revision: Jan 2003, WDH
00036   const int maxnOscillStep = 5 ;
00037   const int maxnDivergingStep = 5 ;
00038 
00039   // initialize
00040   _status = TrkErrCode::succeed;
00041   _flt1 = f1;
00042   _flt2 = f2;
00043   
00044   static HepPoint3D newPos1, newPos2 ;
00045   double delta(0), prevdelta(0) ;
00046   int nOscillStep(0) ;
00047   int nDivergingStep(0) ;
00048   bool finished = false ;
00049   int istep(0) ;
00050 
00051   for (istep=0; istep < _maxTry && !finished; ++istep) {
00052     double prevflt1  = flt1();
00053     double prevflt2  = flt2() ;
00054     double prevprevdelta = prevdelta ;
00055     prevdelta = delta;
00056     
00057     stepTowardPoca(traj1, traj2);
00058     if( status().failure() ) {
00059       // failure in stepTowardPoca
00060       finished=true ;
00061     } else {
00062       newPos1 = traj1.position(flt1());
00063       newPos2 = traj2.position(flt2());
00064       delta = (newPos1 - newPos2).mag();
00065       double step1 = flt1() - prevflt1;
00066       double step2 = flt2() - prevflt2;
00067       int pathDir1 = (step1 > 0.) ? 1 : -1;
00068       int pathDir2 = (step2 > 0.) ? 1 : -1;
00069       
00070       // Can we stop stepping?
00071       double distToErr1 = traj1.distTo1stError(prevflt1, precision(), pathDir1);
00072       double distToErr2 = traj2.distTo1stError(prevflt2, precision(), pathDir2);
00073       
00074       // converged if very small steps, or if parallel
00075       finished = 
00076         (fabs(step1) < distToErr1 && fabs(step2) < distToErr2 ) ||
00077         (status().success() == 3) ;
00078  
00079       // we have to catch some problematic cases
00080       if( !finished && istep>2 && delta > prevdelta) {
00081         if( prevdelta > prevprevdelta) {
00082           // diverging
00083           if(++nDivergingStep>maxnDivergingStep) { 
00084             _status.setFailure(2) ; // code for `Failed to converge'
00085             finished = true ;
00086           }
00087         } else {
00088           nDivergingStep=0;
00089           // oscillating
00090           if(++nOscillStep>maxnOscillStep) {
00091             // bail out of oscillation. since the previous step was
00092             // better, use that one.
00093             _flt1 = prevflt1 ;
00094             _flt2 = prevflt2 ;
00095             _status.setSuccess(21, "Oscillating poca.") ;
00096             finished = true ;
00097           } else {
00098             // we might be oscillating, but we could also just have
00099             // stepped over the minimum. choose a solution `in
00100             // between'.
00101             _flt1 = prevflt1 + 0.5*step1 ;
00102             _flt2 = prevflt2 + 0.5*step2 ;          
00103             newPos1 = traj1.position(_flt1) ;
00104             newPos2 = traj2.position(_flt2) ;
00105             delta = (newPos1 - newPos2).mag() ;
00106           }
00107         }
00108       } 
00109     }
00110   }
00111   if(!finished) _status.setSuccess(2) ; // code for 'not converged' (yet)
00112 }

double TrkPocaBase::precision (   )  [inline, inherited]
 

double TrkPocaBase::precision (   )  [inline, inherited]
  

00059 {return _precision;}

const TrkErrCode& TrkPocaBase::status (   )  const [inline, inherited]
 

const TrkErrCode & TrkPocaBase::status (   )  const [inline, inherited]
  

00062 {return _status;}

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

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

00251 {
00252   // Unsightly uninitialized variables:
00253   static Hep3Vector dir, delDir;
00254   static HepPoint3D trajPos;
00255 
00256   traj.getInfo(flt1(), trajPos, dir, delDir);
00257   Hep3Vector delta = ((CLHEP::Hep3Vector)trajPos) - ((CLHEP::Hep3Vector)pt);
00258   double denom = 1. + delta.dot(delDir);
00259   if (fabs(denom)*_maxDist < 1. ) {
00260     _status.setFailure(11, "TrkPoca::ambiguous tight looper.");
00261     return;
00262   }
00263   double df = -delta.dot(dir) / fabs(denom);
00264   int pathDir = (df > 0.) ? 1 : -1;
00265 
00266   // Don't try going farther than worst parabolic approximation will allow:
00267   double distToErr = traj.distTo2ndError(flt1(), _extrapToler, pathDir);
00268   if (fabs(df)>distToErr) df = (df>0?distToErr:-distToErr);
00269   // Make the step slightly longer -- prevents quitting at kinks
00270   df += 0.001 * pathDir * precision();
00271   _flt1 += df;
00272 }

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

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

00177 { 
00178   // Last revision: Jan 2003, WDH
00179   
00180   // A bunch of unsightly uninitialized variables:
00181   static Hep3Vector dir1, dir2;
00182   static Hep3Vector delDir1, delDir2;
00183   static HepPoint3D pos1, pos2;
00184 
00185   traj1.getInfo(flt1(), pos1, dir1, delDir1);
00186   traj2.getInfo(flt2(), pos2, dir2, delDir2);
00187   Hep3Vector delta = ((CLHEP::Hep3Vector)pos1) - ((CLHEP::Hep3Vector)pos2);
00188   double ua = -delta.dot(dir1);
00189   double ub =  delta.dot(dir2);
00190   double caa = dir1.mag2() + delta.dot(delDir1);
00191   double cbb = dir2.mag2() - delta.dot(delDir2);
00192   double cab = -dir1.dot(dir2);
00193   double det = caa * cbb - cab * cab;
00194 
00195   if(det<0) {
00196     // get rid of second order terms
00197     caa = dir1.mag2() ;
00198     cbb = dir2.mag2() ;
00199     det = caa * cbb - cab * cab;
00200   }
00201           
00202   if ( det < 1.e-8) {
00203     // If they are parallel (in quadratic approximation) give up
00204     _status.setSuccess(3);
00205     return;
00206   }
00207 
00208   double df1 = (ua * cbb - ub * cab)/det;
00209   int pathDir1 = (df1 > 0) ? 1 : -1;
00210   double df2 = (ub * caa - ua * cab)/det;
00211   int pathDir2 = (df2 > 0) ? 1 : -1;
00212 
00213   // Don't try going farther than worst parabolic approximation will
00214   // allow: Since ` _extrapToler' is large, this cut effectively only
00215   // takes care that we don't make large jumps past the kink in a
00216   // piecewise trajectory.
00217 
00218   double distToErr1 = traj1.distTo2ndError(flt1(), _extrapToler, pathDir1);
00219   double distToErr2 = traj2.distTo2ndError(flt2(), _extrapToler, pathDir2);
00220 
00221   // Factor to push just over border of piecewise traj (essential!)
00222   const double smudge = 1.01 ; 
00223   if( fabs(df1) > smudge*distToErr1 ) {
00224     // choose solution for which df1 steps just over border
00225     df1 = smudge*distToErr1 * pathDir1 ;
00226     // now recalculate df2, given df1:
00227     df2 = (ub - df1*cab)/cbb ;
00228   }
00229 
00230   if( fabs(df2) > smudge*distToErr2 ) {
00231     // choose solution for which df2 steps just over border
00232     df2 = smudge*distToErr2 * pathDir2 ;
00233     // now recalculate df1, given df2:
00234     df1 = (ua - df2*cab)/cbb ;
00235     // if still not okay,
00236     if( fabs(df1) > smudge*distToErr1 ) {
00237       df1 = smudge*distToErr1 * pathDir1 ;
00238     }
00239   }
00240 
00241   _flt1 += df1;
00242   _flt2 += df2;
00243 
00244   // another check for parallel trajectories
00245   if (fabs(flt1()) > _maxDist && fabs(flt2()) > _maxDist) 
00246     _status.setSuccess(3) ;
00247 }

Member Data Documentation

DifNumber TrkDifPoca::_doca [private]
 

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

double TrkPocaBase::_flt1 [protected, inherited]
 

double TrkPocaBase::_flt2 [protected, inherited]
 

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

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

double TrkPocaBase::_precision [protected, inherited]
 

TrkErrCode TrkPocaBase::_status [protected, inherited]
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 19:06:42 2011 for BOSS6.5.5 by  doxygen 1.3.9.1