#include <TrkPocaBase.h>
Inheritance diagram for TrkPocaBase:
Public Member Functions | |
const TrkErrCode & | status () const |
double | flt1 () const |
double | flt2 () const |
double | precision () |
Protected Member Functions | |
TrkPocaBase (double flt1, double flt2, double precision) | |
TrkPocaBase (double flt1, double precision) | |
virtual | ~TrkPocaBase () |
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. |
Definition at line 25 of file TrkPocaBase.h.
TrkPocaBase::TrkPocaBase | ( | double | flt1, | |
double | flt2, | |||
double | precision | |||
) | [protected] |
Definition at line 25 of file TrkPocaBase.cxx.
00026 : _precision(prec), _flt1(f1), _flt2(f2),_status(TrkErrCode::fail) 00027 { 00028 assert(prec > 0.); 00029 }
TrkPocaBase::TrkPocaBase | ( | double | flt1, | |
double | precision | |||
) | [protected] |
Definition at line 114 of file TrkPocaBase.cxx.
00115 : _precision(prec), _flt1(f1), _flt2(0), _status(TrkErrCode::fail) 00116 { 00117 }
TrkPocaBase::~TrkPocaBase | ( | ) | [protected, virtual] |
double TrkPocaBase::flt1 | ( | ) | const [inline] |
Definition at line 65 of file TrkPocaBase.h.
References _flt1.
Referenced by TrkDifPieceTraj::append(), TrkPoca::calcDist(), TrkDifPoca::calcDist(), TrkSimpTraj::changePoint(), TrkPocaXY::fltl1(), TrkSimpleRep::getAllWeights(), TrkCompTrk::getAllWeights(), minimize(), TrkDifPieceTraj::prepend(), stepToPointPoca(), stepTowardPoca(), TrkPoca::TrkPoca(), TrkPocaXY::TrkPocaXY(), and TrkHitOnTrk::updatePoca().
00065 {return _flt1;}
double TrkPocaBase::flt2 | ( | ) | const [inline] |
Definition at line 68 of file TrkPocaBase.h.
References _flt2.
Referenced by TrkPoca::calcDist(), TrkDifPoca::calcDist(), MdcUtilitySvc::docaPatPar(), TrkPocaXY::fltl2(), minimize(), stepTowardPoca(), and TrkHitOnTrk::updatePoca().
00068 {return _flt2;}
void TrkPocaBase::minimize | ( | const Trajectory & | traj1, | |
double | f1, | |||
const HepPoint3D & | pt | |||
) | [protected] |
Definition at line 120 of file TrkPocaBase.cxx.
References _flt1, _maxTry, _status, Trajectory::distTo1stError(), flt1(), genRecEmupikp::i, Trajectory::position(), precision(), TrkErrCode::setFailure(), status(), stepToPointPoca(), and TrkErrCode::succeed.
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] |
Definition at line 32 of file TrkPocaBase.cxx.
References _flt1, _flt2, _maxTry, _status, Trajectory::distTo1stError(), false, flt1(), flt2(), Trajectory::position(), precision(), TrkErrCode::setFailure(), TrkErrCode::setSuccess(), status(), stepTowardPoca(), TrkErrCode::succeed, and true.
Referenced by TrkDifPoca::TrkDifPoca(), and TrkPoca::TrkPoca().
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] |
Definition at line 59 of file TrkPocaBase.h.
References _precision.
Referenced by TrkDifPoca::calcDist(), minimize(), and stepToPointPoca().
00059 {return _precision;}
const TrkErrCode & TrkPocaBase::status | ( | ) | const [inline] |
Definition at line 62 of file TrkPocaBase.h.
References _status.
Referenced by TrkDifPieceTraj::append(), TrkPoca::calcDist(), TrkDifPoca::calcDist(), TrkSimpTraj::changePoint(), MdcHitOnTrack::dcaToWire(), TrkHitOnTrk::getFitStuff(), minimize(), TrkDifPieceTraj::prepend(), TrkDifPoca::TrkDifPoca(), TrkPoca::TrkPoca(), TrkPocaXY::TrkPocaXY(), and TrkHitOnTrk::updatePoca().
00062 {return _status;}
void TrkPocaBase::stepToPointPoca | ( | const Trajectory & | traj, | |
const HepPoint3D & | pt | |||
) | [protected] |
Definition at line 250 of file TrkPocaBase.cxx.
References _extrapToler, _flt1, _maxDist, _status, Trajectory::distTo2ndError(), flt1(), Trajectory::getInfo(), precision(), and TrkErrCode::setFailure().
Referenced by minimize().
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] |
Definition at line 176 of file TrkPocaBase.cxx.
References _extrapToler, _flt1, _flt2, _maxDist, _status, Trajectory::distTo2ndError(), flt1(), flt2(), Trajectory::getInfo(), and TrkErrCode::setSuccess().
Referenced by minimize().
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 }
double TrkPocaBase::_extrapToler = 2. [static, protected] |
double TrkPocaBase::_flt1 [protected] |
Definition at line 40 of file TrkPocaBase.h.
Referenced by flt1(), minimize(), stepToPointPoca(), stepTowardPoca(), and TrkPocaXY::TrkPocaXY().
double TrkPocaBase::_flt2 [protected] |
Definition at line 41 of file TrkPocaBase.h.
Referenced by flt2(), minimize(), stepTowardPoca(), and TrkPocaXY::TrkPocaXY().
double TrkPocaBase::_maxDist = 1.e7 [static, protected] |
int TrkPocaBase::_maxTry = 500 [static, protected] |
double TrkPocaBase::_precision [protected] |
TrkErrCode TrkPocaBase::_status [protected] |
Definition at line 42 of file TrkPocaBase.h.
Referenced by TrkDifPoca::calcDist(), TrkPocaXY::interLineCircle(), minimize(), status(), stepToPointPoca(), stepTowardPoca(), and TrkPocaXY::TrkPocaXY().