BOSS 7.0.4
BESIII Offline Software System
Loading...
Searching...
No Matches
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
 
 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 TrkErrCodestatus () const
 
double flt1 () const
 
double flt2 () const
 
double precision ()
 
const TrkErrCodestatus () 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)
 
 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

Constructor & Destructor Documentation

◆ TrkPocaXY() [1/4]

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

Definition at line 23 of file TrkPocaXY.cxx.

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

◆ TrkPocaXY() [2/4]

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.

46 : TrkPocaBase(flt1,flt2,prec) , _docaxy(-9999.0) {
47
48 // this traj-traj version starts by approximating the trejectories as
49 // either a line or a circle and treats separately the line-line,
50 // circle-circle, and line-circle cases.
51
52 _flt1 = flt1;
53 _flt2 = flt2;
54 double delta1(10*prec);
55 double delta2(10*prec);
56 unsigned niter(0);
57 static const unsigned maxiter(20);
58 while(niter < maxiter &&
59 ( delta1 > prec || delta2 > prec ) ){
60 // get positions and directions
61 HepPoint3D pos1 = traj1.position(_flt1);
62 HepPoint3D pos1b;
63 if(niter == 0) pos1b = pos1;
64 HepPoint3D pos2 = traj2.position(_flt2);
65 Hep3Vector dir1 = traj1.direction(_flt1);
66 Hep3Vector dir2 = traj2.direction(_flt2);
67 Hep3Vector dd1 = traj1.delDirect(_flt1);
68 Hep3Vector dd2 = traj2.delDirect(_flt2);
69 double curv1 = traj1.curvature(_flt1);
70 double curv2 = traj2.curvature(_flt2);
71 double m1, m2, q1, q2;
72 double r1, r2, xc1,xc2,yc1,yc2,sr1,sr2;
73 if(curv1 == 0) {
74 // convert to m*x+q
75 if(dir1[0] == 0){
76 m1 = 1.0e12;
77 }else{
78 m1 = dir1[1]/dir1[0];
79 }
80 q1 = pos1.y()-pos1.x()*m1;
81 }else{
82 // get circle parameters
83 r1 = (1- dir1[2]*dir1[2])/curv1;
84 sr1=r1;
85 if(dir1[0]*dd1[1] - dir1[1]*dd1[0] < 0) sr1 = -r1;
86 double cosphi1 = dir1[0]/sqrt(dir1[0]*dir1[0]+dir1[1]*dir1[1]);
87 double sinphi1 = cosphi1 * dir1[1]/dir1[0];
88 xc1 = pos1.x() - sr1 * sinphi1;
89 yc1 = pos1.y() + sr1 * cosphi1;
90 }
91 if(curv2 == 0) {
92 if(dir2[0] == 0){
93 m2 = 1.0e12;
94 }else{
95 m2 = dir2[1]/dir2[0];
96 }
97 q2 = pos2.y()-pos2.x()*m2;
98 }else{
99 r2 = (1-dir2[2]*dir2[2])/curv2;
100 sr2 = r2;
101 if(dir2[0]*dd2[1] - dir2[1]*dd2[0] < 0) sr2 = -r2;
102 double cosphi2 = dir2[0]/sqrt(dir2[0]*dir2[0]+dir2[1]*dir2[1]);
103 double sinphi2 = cosphi2 * dir2[1]/dir2[0];
104 xc2 = pos2.x() - sr2 * sinphi2;
105 yc2 = pos2.y() + sr2 * cosphi2;
106 }
107 double xint, yint, xint1, xint2, yint1, yint2, absdoca;
108 _docaxy = 0;
110
111 // First the line-line case
112
113 if(curv1==0 && curv2==0){
114 // do the intersection in 2d
115 interTwoLines(m1, q1, m2, q2, xint, yint);
116 }
117
118 // next do the two circle case
119
120 if(curv1 !=0 && curv2 !=0){
121 //There are four cases
122 double cdist = sqrt((xc1-xc2)*(xc1-xc2)+(yc1-yc2)*(yc1-yc2));
123 // a - coincident centers
124 if (fabs(cdist) < 1.e-12 ) {
125 //the algorithm fails because the points have all the same distance
127 "TrkPocaXY:: the two circles have the same center...");
128 return;
129 }
130 // b - Actual intersection
131 if ( (fabs(r1-r2) < cdist) && (cdist < r1+r2 ) ) {
132 interTwoCircles(xc1,yc1,r1,xc2,yc2,r2,xint1,yint1,xint2,yint2);
133 double dist1 = (xint1-pos1b.x())*(xint1-pos1b.x()) +
134 (yint1-pos1b.y())*(yint1-pos1b.y());
135 double dist2 = (xint2-pos1b.x())*(xint2-pos1b.x()) +
136 (yint2-pos1b.y())*(yint2-pos1b.y());
137
138 //choose closest to begining of track1
139 if(dist1<dist2){
140 xint = xint1;
141 yint = yint1;
142 } else {
143 xint = xint2;
144 yint = yint2;
145 }
146 }
147
148 // c - nested circles and d - separated circles
149
150 if ( (fabs(r1-r2) > cdist) || // nested circles
151 ( cdist > (r1+r2) )) { // separated circles
152 // line going through the centers of the two circles
153
154 double m = (yc1-yc2)/(xc1-xc2); // y = m*x+q
155 double q = yc1 - m*xc1;
156
157 // intersection points between the line and the two circles
158
159 double xint1, yint1, xint2, yint2, zOfDCrossT1;
160
161 interLineCircle(m, q, xc1, yc1, r1, xint1, yint1, xint2, yint2);
162
163 double xint3, yint3, xint4, yint4 ;
164 interLineCircle(m, q, xc2, yc2, r2, xint3, yint3, xint4, yint4);
165 if (fabs(r1-r2) > cdist) { // nested circles
166 absdoca = fabs(r1-r2)-cdist;
167#ifdef MDCPATREC_DEBUG
168 cout << "ErrMsg(debugging) doing nested circles in TrkPocaXY " << endl;
169#endif
170
171 double dist1_3 = pow((xint1-xint3),2.) + pow((yint1-yint3),2.);
172 double dist2_4 = pow((xint2-xint4),2.) + pow((yint2-yint4),2.);
173
174 if(dist1_3<dist2_4) {
175 xint = 0.5*(xint1+xint3);
176 yint = 0.5*(yint1+yint3);
177 zOfDCrossT1 = (xint3-xint1)*dir1[1]-(yint3-yint1)*dir1[0];
178 } else {
179 xint = 0.5*(xint2+xint4);
180 yint = 0.5*(yint2+yint4);
181 zOfDCrossT1 = (xint4-xint2)*dir1[1]-(yint4-yint2)*dir1[0];
182 }
183 _docaxy = absdoca;
184 if( zOfDCrossT1 > 0) _docaxy = -absdoca;
185 }
186 if( cdist > (r1+r2) ) { //separated circles
187 absdoca = cdist - (r1+r2);
188#ifdef MDCPATREC_DEBUG
189 cout << "ErrMsg(debugging) doing separated circles in TrkPocaXY"<< endl;
190#endif
191
192 double dist2_3 = pow((xint2-xint3),2.) + pow((yint2-yint3),2.);
193 double dist1_4 = pow((xint1-xint4),2.) + pow((yint1-yint4),2.);
194 if (dist2_3<dist1_4){
195 xint = 0.5*(xint2+xint3);
196 yint = 0.5*(yint2+yint3);
197 zOfDCrossT1 = (xint3-xint1)*dir1[1]-(yint3-yint1)*dir1[0];
198 } else {
199 xint = 0.5*(xint1+xint4);
200 yint = 0.5*(yint1+yint4);
201 zOfDCrossT1 = (xint4-xint2)*dir1[1]-(yint4-yint2)*dir1[0];
202 }
203 _docaxy = absdoca;
204 if( zOfDCrossT1 > 0) _docaxy = -absdoca;
205 }
206 }
207 }
208
209 // Now do the line-circle case
210
211 if((curv1 == 0 && curv2 !=0) || (curv1 != 0 && curv2 == 0)) {
212 // distance between the line and the circle center
213 HepVector dirT;
214 double m,q,r,xc,yc, zOfDCrossT1;
215 if(curv1==0) {m=m1; q=q1; r=r2; xc=xc2; yc=yc2; dirT=dir2;}
216 else {m=m2; q=q2; r=r1; xc=xc1; yc=yc1; dirT=dir1;}
217
218 double dist = fabs((m*xc-yc+q)/sqrt(1+m*m));
219 if(dist <= r) {
220
221 // the intersection points
222
223 interLineCircle(m,q,xc,yc,r, xint1, yint1, xint2, yint2);
224 double dist1 = (xint1-pos1b.x())*(xint1-pos1b.x()) +
225 (yint1-pos1b.y())*(yint1-pos1b.y());
226 double dist2 = (xint2-pos1b.x())*(xint2-pos1b.x()) +
227 (yint2-pos1b.y())*(yint2-pos1b.y());
228 //choose closest to the beginning of track 1
229 if(dist1<dist2){
230 xint = xint1;
231 yint = yint1;
232 } else {
233 xint = xint2;
234 yint = yint2;
235 }
236 } else { // no intersection points
237#ifdef MDCPATREC_DEBUG
238 cout << "ErrMsg(debugging) doing separated line-circle in TrkPocaXY"<<endl;
239#endif
240
241 // line going through the circle center and perpendicular to traj1
242
243 double mperp = -1./m;
244 double qperp = yc - mperp*xc;
245
246 // intersection between this line and the two trajectories
247
248 interLineCircle(mperp, qperp, xc, yc, r, xint1,yint1,xint2,yint2);
249
250 double xint3,yint3;
251 interTwoLines(m, q, mperp, qperp, xint3, yint3);
252
253 double dist1_3 = pow((xint1-xint3),2.) + pow((yint1-yint3),2.);
254 double dist2_3 = pow((xint2-xint3),2.) + pow((yint2-yint3),2.);
255 if (dist1_3<dist2_3) {
256 xint = 0.5*(xint3 + xint1);
257 yint = 0.5*(yint3 + yint1);
258 absdoca = sqrt((xint3-xint1)*(xint3-xint1)+
259 (yint3-yint1)*(yint3-yint1));
260 zOfDCrossT1 = (xint3-xint1)*dirT[1]-(yint3-yint1)*dirT[0];
261 }
262 else {
263 xint = 0.5*(xint3 + xint2);
264 yint = 0.5*(yint3 + yint2);
265 absdoca = sqrt((xint3-xint2)*(xint3-xint2)+
266 (yint3-yint2)*(yint3-yint2));
267 zOfDCrossT1 = (xint3-xint2)*dirT[1]-(yint3-yint2)*dirT[0];
268 }
269 _docaxy = absdoca;
270 if( zOfDCrossT1 > 0) _docaxy = -absdoca;
271 }
272 }
273
274 // get the flight lengths for the intersection point
275
276 HepPoint3D intpt( xint, yint, 0.0);
277 TrkPocaXY poca1(traj1,_flt1,intpt,prec);
278 TrkPocaXY poca2(traj2,_flt2,intpt,prec);
279 _status = poca2.status();
280 if(poca1.status().success() && poca2.status().success()) {
281 delta1 = fabs(_flt1 - poca1.flt1());
282 delta2 = fabs(_flt2 - poca2.flt1());
283 _flt1 = poca1.flt1();
284 _flt2 = poca2.flt1();
285 }else{
286#ifdef MDCPATREC_WARNING
287 cout << "ErrMsg(warning)" << " poca fialure " << endl;
288#endif
289 if(poca1.status().failure()) _status=poca1.status();
290 break;
291 }
292 niter++;
293 }
294#ifdef MDCPATREC_DEBUG
295 cout <<"ErrMsg(debugging)" <<"TrkPocaXY iterations = " << niter-1
296 << " flight lengths " << _flt1 <<" " << _flt2 << endl
297 << " deltas " << delta1 <<" " << delta2 << endl;
298#endif
299 if(niter == maxiter){
300#ifdef MDCPATREC_ROUTINE
301 cout << "ErrMsg(routine)" << "TrkPocaXY:: did not converge" << endl;
302#endif
303 _status.setFailure(13, "TrkPocaXY:: did not converge");
304 }
305}
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition: KKsem.h:33
virtual HepPoint3D position(double) const =0
virtual Hep3Vector delDirect(double) const =0
virtual Hep3Vector direction(double) const =0
virtual double curvature(double) const =0
void setSuccess(int i, const char *str=0)
void setFailure(int i, const char *str=0)

◆ ~TrkPocaXY() [1/2]

TrkPocaXY::~TrkPocaXY ( )
inline

Definition at line 43 of file InstallArea/include/TrkBase/TrkBase/TrkPocaXY.h.

43{}

◆ TrkPocaXY() [3/4]

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

◆ TrkPocaXY() [4/4]

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

◆ ~TrkPocaXY() [2/2]

TrkPocaXY::~TrkPocaXY ( )
inline

Member Function Documentation

◆ docaXY() [1/2]

double TrkPocaXY::docaXY ( ) const
inline

Definition at line 73 of file InstallArea/include/TrkBase/TrkBase/TrkPocaXY.h.

73{return _docaxy;}

◆ docaXY() [2/2]

double TrkPocaXY::docaXY ( ) const
inline

◆ fltl1() [1/2]

double TrkPocaXY::fltl1 ( ) const
inline

Definition at line 49 of file InstallArea/include/TrkBase/TrkBase/TrkPocaXY.h.

49{ return flt1(); }

◆ fltl1() [2/2]

double TrkPocaXY::fltl1 ( ) const
inline

◆ fltl2() [1/2]

double TrkPocaXY::fltl2 ( ) const
inline

Definition at line 50 of file InstallArea/include/TrkBase/TrkBase/TrkPocaXY.h.

50{ return flt2(); }

◆ fltl2() [2/2]

double TrkPocaXY::fltl2 ( ) const
inline

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