TrkPocaXY Class Reference

#include <TrkPocaXY.h>

Inheritance diagram for TrkPocaXY:

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
Public Member Functions inherited from TrkPocaBase
const TrkErrCode &	status () const
double	flt1 () const
double	flt2 () const
double	precision ()

Additional Inherited Members
Protected Member Functions inherited from TrkPocaBase
	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 inherited from TrkPocaBase
double	_precision
double	_flt1
double	_flt2
TrkErrCode	_status
Static Protected Attributes inherited from TrkPocaBase
static double	_maxDist = 1.e7
static int	_maxTry = 500
static double	_extrapToler = 2.

Detailed Description

Definition at line 27 of file TrkPocaXY.h.

Constructor & Destructor Documentation

TrkPocaXY() [1/2]

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

Definition at line 23 of file TrkPocaXY.cxx.

  : 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() [2/2]

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.

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

inline

Definition at line 43 of file TrkPocaXY.h.

{}

Member Function Documentation

docaXY()

double TrkPocaXY::docaXY ( ) const

inline

Definition at line 73 of file TrkPocaXY.h.

{return _docaxy;}

fltl1()

double TrkPocaXY::fltl1 ( ) const

inline

Definition at line 49 of file TrkPocaXY.h.

{ return flt1(); }

fltl2()

double TrkPocaXY::fltl2 ( ) const

inline

Definition at line 50 of file TrkPocaXY.h.

{ return flt2(); }

---

The documentation for this class was generated from the following files:

• TrkPocaXY.h
• TrkPocaXY.cxx